#########################
### This code replicates the main tables and plots for Justice as Checks and Balances
### using the dataset for individual calims from the AGN
### Author: Edgar Franco-Vivanco
### Date: May 2021
##########################

### This script contains the code to create
## Figure 1.- Stacked plot for proportion of Documents by Type of Ruling (1590-1820)
## Table  3.- Total and proportions by topics
## Main results: Table 4 and Figure 2 of predicted probabilities
## Merits of the case: Table A11
## Political Context results: Table A9, A10 and A11

# Preliminaries --
rm(list=ls())
gc()

### Libraries
###########################
# Relevant libraries
# Next, we download packages that H2O depends on.
pkgs <- c("plm", "dplyr", "magrittr", "ggplot2", "lmtest", "stargazer",
          "aod", "multiwayvcov", "xtable", "sampleSelection", "car", "emmeans", "phia", "foreign",
          "car", "gridExtra", "RcppRoll"
)
for (pkg in pkgs) {
  if (! (pkg %in% rownames(installed.packages()))) { install.packages(pkg) }
  lapply(pkg, library, character.only=T )
}

# Set path to your own
setwd("C://Users/franc/Dropbox/PhD/Edgar_Tlaxcala/AGN/")

# DONT RUN, Cleaning files
# X "file"              "year_in,year_end"          "topic"             "titulo"            "productor"         "fojas"             "fojas2"   
# "year"              "name"            "decade"  pop_change_d" tierras_d abusos_d  impuestos_d  comunidad_d
# agn_indios <- read.csv('agn_indios3.csv')
# agn_indios <- agn_indios %>% dplyr::select(X, file, NewID, year_in,year_end,topic,titulo,productor,fojas2,year,name,region, decade, pob, pob_ln, pop_change_d, pop_ratio,
#                                     tierras_d ,abusos_d,  impuestos_d  ,comunidad_d, protect_new2, per_encomienda, power)
# 
# ## Create pop decrease dummy
# agn_indios <- agn_indios %>% mutate(pop_decrease=ifelse(pop_change_d==1,0,1),
#                                     pop_decrease=ifelse(is.na(pop_change_d), NA, pop_decrease)
#                                     )
# 
# ###### Identify cases with only one topic
# agn_indios <- agn_indios %>% mutate(tierras_d3 = ifelse(tierras_d==1 & (abusos_d+impuestos_d+comunidad_d==0),1,0),
#                                             abusos_d3  = ifelse(abusos_d==1  & (tierras_d+impuestos_d+comunidad_d==0),1,0),
#                                             impuestos_d3  = ifelse(impuestos_d==1  & (tierras_d+abusos_d+comunidad_d==0),1,0),
#                                             comunidad_d3  = ifelse(comunidad_d==1  & (tierras_d+abusos_d+impuestos_d==0),1,0),
#                                             ##
#                                             theme = ifelse(tierras_d3==1, "Tierras", NA),
#                                             theme = ifelse(abusos_d3==1,   "Mistreatment", theme),
#                                             theme = ifelse(impuestos_d3==1, "Taxes", theme),
#                                             theme = ifelse(comunidad_d3==1, "Community", theme),
#                                             theme = as.factor(theme),
#                                             theme = relevel(theme, ref="Tierras"),
#                                             ##
#                                             only_one = ifelse(abusos_d+impuestos_d+comunidad_d+tierras_d==1,1,0)
# )
# 
# ## Save final
# write.csv(agn_indios, "C://Users/franc/Dropbox/PhD/Edgar_Tlaxcala/AGN/agn_indios_2021.csv")
### READ file
agn_indios <- read.csv("C://Users/franc/Dropbox/PhD/Edgar_Tlaxcala/AGN/agn_indios_2021.csv")

## For the rest of the analysis we need to delete the provinces of Valladolid and Guachinango
## Because the location is ambigous (Two Valladolids and two Guachinangos)
## Also, exclude Mexico City because this was the center of the GIC and many cases not originated there had Mexico City as province
#agn_indios$X[duplicated(agn_indios$X)]
agn_indios_sub <- agn_indios %>% filter(name!="Valladolid", name!="Guachinango", name!="Guanchinango", name!="Mexico") %>% filter(year>1592) 

### Read rebellions data 
rebel <- read.csv("C://Users/franc/Dropbox/PhD/dissertation/strategies_resistance_accomodation/rebel5.csv")
rebel <- rebel %>% dplyr::select(-X) %>% dplyr::select(name, decade, small_scale2,rebel_hist2) 
#write final 
rebel <- write.csv(rebel, "C://Users/franc/Dropbox/PhD/dissertation/strategies_resistance_accomodation/rebel_2021.csv")
rebel <- read.csv("C://Users/franc/Dropbox/PhD/dissertation/strategies_resistance_accomodation/rebel_2021.csv")

### Merge shorter dataset with panel 
agn_indios_sub <- agn_indios_sub  %>% left_join(rebel, by=c("name", "decade"))

##### Distance to Mexico city
distance_mx <- read.dbf("C://Users/franc/Dropbox/Aztec/Gerhard/near3.dbf")
mex_ter     <- read.dbf("C://Users/franc/Dropbox/Aztec/Gerhard/gerhard_centroid.dbf")

mex_ter <- mex_ter %>% dplyr::select(ORIG_FID, NewID) %>% rename(IN_FID=ORIG_FID)
distance_mx <- distance_mx %>% dplyr::select(IN_FID, NEAR_DIST)
mex_ter <- mex_ter %>% left_join(distance_mx, by="IN_FID") %>% dplyr::select(-IN_FID)
# Final merge 
agn_indios_sub <- agn_indios_sub  %>% left_join(mex_ter, by=c("NewID")) %>% mutate(distance=NEAR_DIST/100000)


####################
##  
##   Figure 1 . Stacked Plot
##
##################
##Inspired by http://www.r-graph-gallery.com/136-stacked-area-chart/ 

## Group by category
agn_noProtect <- agn_indios %>% filter(protect_new2==0) %>%  group_by( year ) %>% summarise(count=n()) %>% 
  mutate(Topic="Unfavorable ruling")
agn_Protect <- agn_indios %>% filter(protect_new2==1)  %>%  group_by( year ) %>% summarise(count=n()) %>% 
  mutate(Topic="Favorable ruling")

agn_group2        <- as.data.frame(rbind(agn_Protect , agn_noProtect))


## Calculate proportion
my_fun           <- function(vec){ as.numeric(vec[2]) / sum(agn_group2$count[agn_group2$year==vec[1]], na.rm=T) *100 }
agn_group2$prop  <- apply(agn_group2 , 1 , my_fun)
agn_group2$Topic <- as.factor(agn_group2$Topic)
agn_group2$Topic <- factor(agn_group2$Topic, levels(agn_group2$Topic)[c(2,1)])

agn_group2 <- agn_group2 %>% filter(year<1821 & year>1590) %>% rename(Resolution=Topic)

agn_period_plot2 <- ggplot(agn_group2, aes(year, prop, fill=Resolution)) + geom_area() + scale_fill_manual(values=c("gray45", "gray70"))
agn_period_plot2 <- agn_period_plot2 + theme_bw() + xlab("Year") + ylab("Percentage of Documents")
agn_period_plot2 <- agn_period_plot2 +  theme(legend.position="bottom", legend.title = element_blank(),
                                              legend.box.background = element_rect(colour = "black"),
                                              panel.border = element_rect(colour = "black", fill=NA)
                                              )
agn_period_plot2

ggsave(agn_period_plot2, file="C:/Users/franc/Dropbox/PhD/dissertation/paper_colonialism_liberalims/tex/FINAL_version_2021/figures_paper/figure1.tiff", width=8)

## Calculate averages to mention in text
mean(agn_indios$protect_new2, na.rm=T) # Percentage protection
mean(agn_indios$protect_new2[which(agn_indios$year>=1600 & agn_indios$year<1700)]) #17th century
mean(agn_indios$protect_new2[which(agn_indios$year>=1700 & agn_indios$year<1800)]) # 18th

########################################################################
##  
##   Table 3: Totals by topic and Proportions
##
########################################################################

### Land Conflicts
table(agn_indios$tierras_d)[2]
table(agn_indios$tierras_d, agn_indios$protect_new2)/ table(agn_indios$tierras_d)[2]

## Mistreatment
table(agn_indios$abusos_d)[2]
table(agn_indios$abusos_d, agn_indios$protect_new2)/ table(agn_indios$abusos_d)[2]

## Taxes
table(agn_indios$impuestos_d)[2]
table(agn_indios$impuestos_d, agn_indios$protect_new2)/ table(agn_indios$impuestos_d)[2]

## Community
table(agn_indios$comunidad_d)[2]
table(agn_indios$comunidad_d, agn_indios$protect_new2)/ table(agn_indios$comunidad_d)[2]

## Protection
table(agn_indios$protect_new2)[2]/sum(table(agn_indios$protect_new2))

#################################################################
##
##   Main results Topics and Balance of Power
##   Table 4 (main text)
##
#################################################################

#### Table 4 --------------------
## Province and Fixed effecs #####################
#### All sample (unconditional)
model_dict_all <- lm(protect_new2 ~tierras_d+ abusos_d +impuestos_d + comunidad_d+as.factor(year)+name*year,
                     data=agn_indios_sub)
model_dict_all_vcov <- cluster.vcov(model_dict_all   , ~ agn_indios_sub$name)
model_dict_all_vcov <- coeftest(model_dict_all , model_dict_all_vcov)


#### Conditioned on Population Decrease
model_dict_pop_dec <- lm(protect_new2 ~tierras_d+ abusos_d +impuestos_d + comunidad_d+ as.factor(year)+name*year,
                         data=agn_indios_sub[which(agn_indios_sub$pop_ratio<1),])
model_dict_pop_dec_vcov <- cluster.vcov(model_dict_pop_dec  , ~ agn_indios_sub[which(agn_indios_sub$pop_ratio<1),]$name)
model_dict_pop_dec_vcov <- coeftest(model_dict_pop_dec , model_dict_pop_dec_vcov)

#### Population Increase 
model_dict_pop_inc <- lm(protect_new2 ~tierras_d+ abusos_d +impuestos_d + comunidad_d + as.factor(year)+name*year,
                         data=agn_indios_sub[which(agn_indios_sub$pop_ratio>=1),])

model_dict_pop_inc_vcov <- cluster.vcov(model_dict_pop_inc  , ~ agn_indios_sub[which(agn_indios_sub$pop_ratio>=1),]$name)
model_dict_pop_inc_vcov <- coeftest(model_dict_pop_inc , model_dict_pop_inc_vcov)


#### LE Power#################
#### LE Strong
## Poulation in Encomienda mean 
enc_mean <- mean(agn_indios_sub$per_encomienda, na.rm=T) # 0.1529982

model_dict_le_strong <- lm(protect_new2 ~tierras_d+ abusos_d +impuestos_d +comunidad_d+as.factor(year)+name*year,
                           data=agn_indios_sub[which(agn_indios_sub$per_encomienda>=enc_mean),])

model_dict_le_strong_vcov <- cluster.vcov(model_dict_le_strong , ~ agn_indios_sub[which(agn_indios_sub$per_encomienda>=enc_mean),]$name)
model_dict_le_strong_vcov <- coeftest(model_dict_le_strong , model_dict_le_strong_vcov)


### LE weak
model_dict_le_weak <- lm(protect_new2 ~tierras_d+ abusos_d +impuestos_d +comunidad_d+as.factor(year)+name*year,
                         data=agn_indios_sub[which(agn_indios_sub$per_encomienda<enc_mean 
                         ),])
model_dict_le_weak_vcov <- cluster.vcov(model_dict_le_weak , ~ agn_indios_sub[which(agn_indios_sub$per_encomienda<enc_mean),]$name)
model_dict_le_weak_vcov <- coeftest(model_dict_le_weak  , model_dict_le_weak_vcov)


#### Conditioned by pop Decrease AND...
## ...LE strong
model_dict_pop_dec_strong <- lm(protect_new2 ~tierras_d+ abusos_d +impuestos_d +comunidad_d+as.factor(year)+name*year,
                                data=agn_indios_sub[which(agn_indios_sub$pop_ratio<1 & agn_indios_sub$per_encomienda>=enc_mean
                                ),])

model_dict_pop_dec_strong_vcov <- cluster.vcov(model_dict_pop_dec_strong , ~ agn_indios_sub[which(agn_indios_sub$pop_ratio<1 & agn_indios_sub$per_encomienda>=enc_mean),]$name)
model_dict_pop_dec_strong_vcov <- coeftest(model_dict_pop_dec_strong , model_dict_pop_dec_strong_vcov)


## ...LE weak
model_dict_pop_dec_weak <- lm(protect_new2 ~tierras_d+ abusos_d +impuestos_d +comunidad_d+as.factor(year)+name*year,
                              data=agn_indios_sub[which( agn_indios_sub$pop_ratio<1 & agn_indios_sub$per_encomienda<enc_mean
                              ),])

model_dict_pop_dec_weak_vcov  <- cluster.vcov(model_dict_pop_dec_weak  , ~ agn_indios_sub[which(agn_indios_sub$pop_ratio<1 & agn_indios_sub$per_encomienda<enc_mean),]$name)
model_dict_pop_dec_weak_vcov<- coeftest(model_dict_pop_dec_weak  ,model_dict_pop_dec_weak_vcov)


#### Tables (Main)
setwd("C://Users/franc/Dropbox/PhD/dissertation/paper_colonialism_liberalims/FINAL_version_2021/new_version/final_tables")

##### Number of Units by regression
all   <- length(unique(agn_indios_sub$name))

dec  <- length(unique(agn_indios_sub$name[which(agn_indios_sub$pop_ratio<1 )]))

inc  <- length(unique(agn_indios_sub$name[which(agn_indios_sub$pop_ratio>=1 )]))

strg  <- length(unique(agn_indios_sub$name[which(agn_indios_sub$per_encomienda>=enc_mean )]))

weak <- length(unique(agn_indios_sub$name[which(agn_indios_sub$per_encomienda<enc_mean )]))

dec_str <- length(unique(agn_indios_sub$name[which(agn_indios_sub$pop_ratio<1 & agn_indios_sub$per_encomienda>=enc_mean )]))

dec_weak <- length(unique(agn_indios_sub$name[which(agn_indios_sub$pop_ratio<1 & agn_indios_sub$per_encomienda<enc_mean)]))

unit_line <- c("Units", all, dec, inc, strg, weak, dec_str, dec_weak)

#### Table 4 ############################
stargazer(model_dict_all , model_dict_pop_dec, model_dict_pop_inc, model_dict_le_strong, model_dict_le_weak,
          model_dict_pop_dec_strong , model_dict_pop_dec_weak, 
          se   = list(model_dict_all_vcov[, 2], model_dict_pop_dec_vcov[, 2], 
                      model_dict_pop_inc_vcov[, 2], model_dict_le_strong_vcov[, 2],model_dict_le_weak_vcov[,2] ,
                      model_dict_pop_dec_strong_vcov[,2], model_dict_pop_dec_weak_vcov[,2]
          ) ,
          p    = list(model_dict_all_vcov[, 4], model_dict_pop_dec_vcov[, 4], 
                      model_dict_pop_inc_vcov[, 4], model_dict_le_strong_vcov[, 4],model_dict_le_weak_vcov[,4] ,
                      model_dict_pop_dec_strong_vcov[,4], model_dict_pop_dec_weak_vcov[,4]
          ),
          dep.var.labels= c("Favorable Ruling"),
          covariate.labels =c("Land Conflicts","Mistreatment",  "Taxes", 
                              "Community"
          ),
          column.labels = c("All", "Population", "Local Elites", "Population Decrease"
          ),
          column.separate = c(1,2,2,2),
          keep = c("tierras_d", "abusos_d", "impuestos_d" , "comunidad_d"), 
          keep.stat = c("n", "adj.rsq"),
          add.lines=list(c("Province FE", "Y", "Y", "Y", "Y", "Y", "Y","Y"),
                         c("Time FE", "Y", "Y", "Y", "Y","Y", "Y","Y") ,
                         c("Province time trends", "Y", "Y", "Y", "Y","Y", "Y","Y") ,
                         unit_line
          ),
          label="topics_balance",
          title =("Favorable Ruling, Topics, and Balance of Powers")
          #out="model_vulner3.tex"
          )

######################################################################
##
##   Figure 2: SIMULATIONS for Predicted Probabilities
##   Note: Zelig is not longer supported in CRAN
#    https://cran.r-project.org/web/packages/Zelig/index.html
##    If you dont have the program installed in your computer or have an
##   updated version of R you can install Zelig following these steps:
##  install.packages("remotes")
##    library(remotes)
##   install_github("cran/Zelig")
library(Zelig)
######################################################################

###### POPULATION
### Pop decrease
z1_1 <- zelig(protect_new2~ tierras_d + abusos_d + impuestos_d + comunidad_d + as.factor(name) +year,
              data=agn_indios_sub[which(agn_indios_sub$pop_ratio<1),] ,  model="logit")

### Pop increase
z1_2 <- zelig(protect_new2~ tierras_d + abusos_d + impuestos_d + comunidad_d+  as.factor(name) + year,
              data=agn_indios_sub[which(agn_indios_sub$pop_ratio>=1 ),], model="logit")

#### LE STRENGHT
### LE strong
# Run line 173 to get enc_mean
z1_3 <- zelig(protect_new2~ tierras_d + abusos_d + impuestos_d + comunidad_d + as.factor(name) + year,
              data=agn_indios_sub[which(agn_indios_sub$per_encomienda>=enc_mean),], model="logit")

### LE weak 
z1_4 <- zelig(protect_new2~ tierras_d + abusos_d + impuestos_d + comunidad_d+  as.factor(name) + year,
              data=agn_indios_sub[which(agn_indios_sub$per_encomienda<enc_mean 
              ),], model="logit")


##### BALANCE of POWERS
## Pop decrease AND local elites strong
z1_5 <- zelig(protect_new2~ tierras_d + abusos_d + impuestos_d + comunidad_d+  as.factor(name) + year,
              data=agn_indios_sub[which(agn_indios_sub$pop_ratio<1 & agn_indios_sub$per_encomienda>=enc_mean
              ),], model="logit")

## Pop decrease and local weak
z1_6 <- zelig(protect_new2~ tierras_d + abusos_d + impuestos_d + comunidad_d+  as.factor(name) + year,
              data=agn_indios_sub[which(agn_indios_sub$pop_ratio<1 & agn_indios_sub$per_encomienda<enc_mean
              ),], model="logit")


#----------------------------------------------#
# Loop to create datasets 
df_dec_pop      <- data.frame()
df_in_pop       <- data.frame()
df_le_str       <- data.frame()
df_le_weak      <- data.frame()
df_dec_pop_le_s <- data.frame()
df_dec_pop_le_w <- data.frame()

v <- list(c(1,0,0,0), c(0,1,0,0), c(0,0,1,0), c(0,0,0,1))
topics <- c("Land", "Mistreatment", "Taxes","Community")

### Loop across topics
for(t in 1:length(topics)){
  ### Decreasing pob
  z1_temp   <- setx(z1_1,tierras_d=v[[t]][1], abusos_d=v[[t]][2], 
                    impuestos_d=v[[t]][3], comunidad_d=v[[t]][4])  
  z1_sim_temp <-  sim(z1_1, x= z1_temp   ) 
  z1_df_temp   <- zelig_qi_to_df(z1_sim_temp) %>% mutate(setx_value=topics[t])
  
  df_dec_pop  <- df_dec_pop  %>% rbind.data.frame(z1_df_temp  %>%  group_by(setx_value) %>% summarise(ev=mean(expected_value), 
                                                                                                      sd=sd(expected_value)))
  
  ### Increasing pob
  z2_temp   <- setx(z1_2, tierras_d=v[[t]][1], abusos_d=v[[t]][2], 
                    impuestos_d=v[[t]][3], comunidad_d=v[[t]][4])  
  z2_sim_temp <-  sim(z1_2, x= z2_temp   ) 
  z2_df_temp   <- zelig_qi_to_df(z2_sim_temp) %>% mutate(setx_value=topics[t])
  
  df_in_pop <- df_in_pop %>% rbind.data.frame(z2_df_temp  %>%  group_by(setx_value) %>% summarise(ev=mean(expected_value), 
                                                                                                  sd=sd(expected_value)))
  
  #### LE strong 
  z3_temp   <- setx(z1_3,  tierras_d=v[[t]][1], abusos_d=v[[t]][2], 
                    impuestos_d=v[[t]][3], comunidad_d=v[[t]][4])  
  z3_sim_temp <-  sim(z1_3, x= z3_temp   ) 
  z3_df_temp   <- zelig_qi_to_df(z3_sim_temp) %>% mutate(setx_value=topics[t])
  
  df_le_str <- df_le_str %>% rbind.data.frame(z3_df_temp  %>%  group_by(setx_value) %>% 
                                                summarise(ev=mean(expected_value), sd=sd(expected_value)))
  
  ### LE weak 
  z4_temp   <- setx(z1_4, tierras_d=v[[t]][1], abusos_d=v[[t]][2], 
                    impuestos_d=v[[t]][3], comunidad_d=v[[t]][4])  
  z4_sim_temp <-  sim(z1_4, x= z4_temp   ) 
  z4_df_temp   <- zelig_qi_to_df(z4_sim_temp) %>% mutate(setx_value=topics[t])
  
  df_le_weak  <- df_le_weak %>% rbind.data.frame(z4_df_temp  %>%  group_by(setx_value) %>% 
                                                   summarise(ev=mean(expected_value), sd=sd(expected_value)))
  
  ### Land pop decline AND local elites strong
  z5_temp   <- setx(z1_5,  tierras_d=v[[t]][1], abusos_d=v[[t]][2], 
                    impuestos_d=v[[t]][3], comunidad_d=v[[t]][4])  
  z5_sim_temp <-  sim(z1_5, x= z5_temp   ) 
  z5_df_temp   <- zelig_qi_to_df(z5_sim_temp) %>% mutate(setx_value=topics[t])
  
  df_dec_pop_le_s <-  df_dec_pop_le_s %>% rbind.data.frame(z5_df_temp  %>%  group_by(setx_value) %>% 
                                                             summarise(ev=mean(expected_value), sd=sd(expected_value)))
  
  ### Land pop decline AND local elites weak
  z6_temp   <- setx(z1_6,  tierras_d=v[[t]][1], abusos_d=v[[t]][2], 
                    impuestos_d=v[[t]][3], comunidad_d=v[[t]][4])  
  z6_sim_temp <-  sim(z1_6, x= z6_temp   ) 
  z6_df_temp   <- zelig_qi_to_df(z6_sim_temp) %>% mutate(setx_value=topics[t])
  
  df_dec_pop_le_w <- df_dec_pop_le_w  %>%  rbind.data.frame(z6_df_temp  %>% group_by(setx_value) %>% 
                                                              summarise(ev=mean(expected_value), sd=sd(expected_value)))
}

#### Join
df_in_pop$pop <- "Increase"; df_dec_pop $pop <- "Decline"
df_le_str$le <- "Strong"; df_le_weak$le <- "Weak"
df_dec_pop_le_s$le <- "Strong"; df_dec_pop_le_w$le <- "Weak"

### df_pop
df_pop  <- rbind.data.frame(df_in_pop,df_dec_pop )
df_pop$id <- rep(c(1,2,3,4),2)
df_pop <- df_pop %>% mutate(pop=as.factor(pop))

### df_le
df_le <- rbind.data.frame(df_le_str,df_le_weak)
df_le$id <- rep(c(1,2,3,4),2)
df_le <- df_le %>% mutate(pop=as.factor(le))

### df pop AND le
df_p_le <- rbind.data.frame(df_dec_pop_le_s,df_dec_pop_le_w)
df_p_le$id <- rep(c(1,2,3,4),2)
df_p_le <- df_p_le  %>% mutate(pop=as.factor(le))

####### Proportions (for text page 35)
#### Change population 
df_pop[1,]
df_pop[5,]

#### Change LE power
df_le[1,]
df_le[5,]

### PLOTS
### GGPLOT Decrease vs Increase
agn_indios_sub_topic <- agn_indios_sub %>% dplyr::select(name, year, protect_new2,pop_ratio, per_encomienda, tierras_d, abusos_d, impuestos_d, comunidad_d) %>% 
  filter((tierras_d==1 | abusos_d==1 | impuestos_d==1 |comunidad_d==1)  & year>=1592,
         agn_indios_sub$name!="Mexico")

rate_dec <- mean(agn_indios_sub_topic$protect_new2[which(agn_indios_sub_topic $pop_ratio<1)], na.rm=T)
rate_inc <- mean(agn_indios_sub_topic$protect_new2[which(agn_indios_sub_topic $pop_ratio>=1)], na.rm=T)

h <- ggplot(df_pop, aes(x=reorder(setx_value, id), y= ev, fill=pop)) + geom_bar(stat="identity",width=0.6, position = "dodge" )
h <- h + geom_pointrange(aes(y= ev, ymin=ev-1.96*sd, ymax=ev+1.96*sd), position = position_dodge(width = 1/2))
h <- h + theme_bw( ) + xlab("")+ ylab("Predicted Probability")+theme(legend.position="bottom", legend.title = element_blank(),
                                                                     legend.box.background = element_rect(colour = "black"),
                                                                     panel.border = element_rect(colour = "black", fill=NA))
h <- h + scale_fill_manual("Population", values = c("Decline" = "grey47", "Increase" = "grey89")) 
h <- h + geom_hline(yintercept = rate_dec, lty="dashed"  ) + geom_hline(yintercept = rate_inc, lty="dotted"  ) 
h <- h + labs(caption="(a) Population Change") +  theme(plot.caption = element_text(hjust=0.5, size=rel(1.2)))
h
ggsave("C:/Users/franc/Dropbox/PhD/dissertation/paper_colonialism_liberalims/tex/FINAL_version_2021/figures_paper/figure2_a.tiff",h, width=6)

### LE power
rate_str  <- mean(agn_indios_sub_topic$protect_new2[which(agn_indios_sub_topic$per_encomienda>=enc_mean)], na.rm=T)
rate_weak <- mean(agn_indios_sub_topic$protect_new2[which(agn_indios_sub_topic$per_encomienda<enc_mean)], na.rm=T)

g <- ggplot(df_le, aes(x=reorder(setx_value, id), y= ev, fill=pop)) + geom_bar(stat="identity",width=0.6, position = "dodge" )
g <- g + geom_pointrange(aes(y= ev, ymin=ev-1.96*sd, ymax=ev+1.96*sd), position = position_dodge(width = 1/2))
g <- g + theme_bw( ) + xlab("")+ ylab("Predicted Probability")+theme(legend.position="bottom", legend.title = element_blank(),
                                                                     legend.box.background = element_rect(colour = "black"),
                                                                     panel.border = element_rect(colour = "black", fill=NA))
g <- g + scale_fill_manual("LE power", values = c("Strong" = "grey47", "Weak" = "grey89"))
g <- g + geom_hline(yintercept = rate_str, lty="dashed"  ) + geom_hline(yintercept = rate_weak, lty="dotted"  ) 
g <- g + labs(caption="(b) Local Elite Power") +  theme(plot.caption = element_text(hjust=0.5, size=rel(1.2)))
g
ggsave("C:/Users/franc/Dropbox/PhD/dissertation/paper_colonialism_liberalims/tex/FINAL_version_2021/figures_paper/figure2_b.tiff",g, width=6)



#############################################
#### Interactions 
#### Table A2.1 in Appendix A2
#### Interactions with pop_change
###########################################

model_dict_int_pop <- lm(protect_new2 ~tierras_d*pop_decrease+ abusos_d*pop_decrease+impuestos_d*pop_decrease + comunidad_d*pop_decrease+
                          as.factor(decade)+name,
                         data=agn_indios_sub)
model_dict_int_pop_vcov <- cluster.vcov(model_dict_int_pop   , ~ agn_indios_sub$name)
model_dict_int_pop_coef <- coeftest(model_dict_int_pop , model_dict_int_pop_vcov)


#######################
#### Interactions with local elites
### All years
model_dict_int_le <- lm(protect_new2 ~tierras_d*power+ abusos_d*power+impuestos_d*power + comunidad_d*power+as.factor(decade)+name,
                        data=agn_indios_sub)
model_dict_int_le_vcov <- cluster.vcov(model_dict_int_le  , ~ agn_indios_sub$name)
model_dict_int_le_coef <- coeftest(model_dict_int_le , model_dict_int_le_vcov )


#####################
#### Triple interaction
#### All
model_dict_int_triple <- lm(protect_new2 ~tierras_d*pop_decrease*power+ abusos_d*pop_decrease*power+impuestos_d*pop_decrease*power + 
                              comunidad_d*pop_decrease*power+as.factor(decade)+name,
                            data=agn_indios_sub)
model_dict_int_triple_vcov <- cluster.vcov(model_dict_int_triple  , ~ agn_indios_sub$name)
model_dict_int_triple_coef <- coeftest(model_dict_int_triple , model_dict_int_triple_vcov )


#### Tables ####
stargazer(model_dict_int_pop, model_dict_int_le, model_dict_int_triple,
          se   = list(model_dict_int_pop_coef[, 2], model_dict_int_le_coef[, 2], 
                      model_dict_int_triple_coef[, 2]
          ) ,
          p    = list(model_dict_int_pop_coef[, 4], model_dict_int_le_coef[, 4], 
                      model_dict_int_triple_coef[, 4]
          ),
          
          ###
          dep.var.labels= c("Favorable Ruling"),
          covariate.labels =c("Land Conflicts", "Land Conflicts x Pop. Decline", "Land Conflicts x Strong LE", "Land Conflicts x Pop. Decline x Strong LE",
                              "Mistreatment", "Mistreatment x Pop. Decline", "Mistreatment x Strong LE", "Mistreatment x Pop. Decline x Strong LE",
                              "Taxes", "Taxes x Pop. Decline" , "Taxes x Strong LE", "Taxes x Pop. Decline x Strong LE",
                              "Community", "Community x Pop. Decline", "Community x Strong LE", "Community x Pop. Decline x Strong LE",
                              "Pop. Decline", "Strong LE" , "Pop. Decline x Strong LE"
          ),
          # column.labels = c("All", "Pre 1700", "Post 1700"
          # ),
          keep = c("tierras_d", "abusos_d", "impuestos_d" , "comunidad_d",
                   "tierras_d:pop_decrease", "pop_decrease:abusos_d", "pop_decrease:impuestos_d",
                   "pop_decrease:comunidad_d", "pop_decrease", 
                   "tierras_d:power", "power:abusos_d", "power:impuestos_d",  "power:comunidad_d", "power"
                   
          ), 
          ##order
          order = c("^tierras_d$", "^tierras_d:pop_decrease$", "^tierras_d:power$", "^tierras_d:pop_decrease:power$",
                    "^abusos_d$","^pop_decrease:abusos_d$", "^power:abusos_d$", "^pop_decrease:power:abusos_d$",
                    "^impuestos_d$", "^pop_decrease:impuestos_d$", "^power:impuestos_d$", "^pop_decrease:power:impuestos_d$", 
                    "^comunidad_d$", "^pop_decrease:comunidad_d$", "^power:comunidad_d$","^pop_decrease:power:comunidad_d$",
                    "^pop_decrease$", "^power$", "^pop_decrease:power$"),
          ##
          keep.stat = c("n", "adj.rsq"),
          add.lines=list(c("Province FE", "Y", "Y", "Y"),
                         c("Time FE", "Y", "Y", "Y") 
          ),
          label="topics_triple_interaction",
          title =("Favorable Ruling, Topics, and Balance of Powers (interactive models)")
          #out="model_vulner3.tex"
)



#########################################
#### Table A2.2
#### Loop through models and topics to collect the population, local elite and double effect
#### Calculate the wald test for linerar equality of coefficients
########################################

model_list <- list(model_dict_int_pop ,model_dict_int_le, model_dict_int_triple )
coef_list  <- list(model_dict_int_pop_coef ,model_dict_int_le_coef, model_dict_int_triple_coef  )
vcov_list  <- list(model_dict_int_pop_vcov ,model_dict_int_le_vcov, model_dict_int_triple_vcov  )
var_list   <- c("tierras_d", "abusos_d", "impuestos_d", "comunidad_d")
  
##
pop_effect    <- rep(NA,12)
power_effect  <- rep(NA,12)
double_effect <- rep(NA,12)
##
p_values_pop    <- rep(NA,12)
p_values_power  <- rep(NA,12)
p_values_double <- rep(NA,12)

#### Begin loop over models
for(m in 1:3){
  #####Begin loop over vars
  for(v in 1:length(var_list)){
  ### Paste to create local variables
  if(var_list[v]=="tierras_d"){
  ### paste topic + pop decrease
  topic_pop   <- paste(var_list[v],":pop_decrease", sep="")
  ### paste topics + power
  topic_pow   <- paste(var_list[v],":power", sep="")
  ### paste double   
  topic_double <-  paste(var_list[v],":pop_decrease:power", sep="")
  
  }else{
    topic_pop   <- paste("pop_decrease:", var_list[v], sep="")
    ### paste topics + power
    topic_pow   <- paste("power:" ,var_list[v], sep="")
    ### paste double   
    topic_double <- paste("pop_decrease:power:",var_list[v], sep="") 
  }
  
  ### create text test
  # test1
  test1 <- paste(var_list[v],"=",var_list[v],"+",topic_pop, "+pop_decrease"  , sep="")
  #test2
  test2 <- paste(var_list[v],"=",var_list[v],"+",topic_pow, "+power"  , sep="")
  #test3 
  test3 <- paste(var_list[v],"=",var_list[v],"+",topic_pop, "+pop_decrease","+",topic_pow, "+power+pop_decrease:power+",topic_double  , sep="")
  
  #### Test significant different (Wald test)
 # https://stats.idre.ucla.edu/r/seminars/interactions-r/#s5a
  
  #### Select in models
  #### Create nv to populate vectors
  if(m==1){
    b1   <- coef_list[[m]][rownames(coef_list[[m]])=="tierras_d",1]
    b2   <- coef_list[[m]][rownames(coef_list[[m]])=="pop_decrease",1]
    b3   <- 0
    b12  <- coef_list[[m]][rownames(coef_list[[m]])==topic_pop,1]
    b13  <- 0
    b123 <- 0
    b23  <- 0
    ##
    nv =v
  }
  if(m==2){
    b1   <- coef_list[[m]][rownames(coef_list[[m]])=="tierras_d",1]
    b2   <- 0
    b3   <- coef_list[[m]][rownames(coef_list[[m]])=="power",1]
    b12  <- 0
    b13  <- coef_list[[m]][rownames(coef_list[[m]])==topic_pow,1]
    b123 <- 0
    b23  <- 0
    ##
    nv=v+4
  }
  if(m==3){
    b1   <- coef_list[[m]][rownames(coef_list[[m]])=="tierras_d",1]
    b2   <- coef_list[[m]][rownames(coef_list[[m]])=="pop_decrease",1]
    b3   <- coef_list[[m]][rownames(coef_list[[m]])=="power",1]
    b12  <- coef_list[[m]][rownames(coef_list[[m]])==topic_pop,1]
    b13  <- coef_list[[m]][rownames(coef_list[[m]])==topic_pow,1]
    b123 <- coef_list[[m]][rownames(coef_list[[m]])==topic_double,1]
    b23  <- coef_list[[m]][rownames(coef_list[[m]])=="pop_decrease:power",1]
    ##
    nv=v+8
    # Double effect
    double_effect[nv] <- (b1+b2+b3+b12+b13+b123+b23)-b1
  }
  ### pop efect
  pop_effect[nv]    <- (b1+b2+b12)-b1
  # power effect
  power_effect[nv]  <- (b1+b3+b13)-b1
  
  
  ##### Collect pvalues from F Wald test
  if(m==1){
    ### Test if land effect alone is different than the effect of land when pop decrease
    p_values_pop[nv]   <- linearHypothesis(model_list[[m]], test1 , vcov=vcov_list[[m]])$`Pr(>F)`[2]
  }
  if(m==2){
    ### Test if land effect alone is different than the effect of land when le strong
    p_values_power[nv]  <- linearHypothesis(model_list[[m]], test2 , vcov=vcov_list[[m]])$`Pr(>F)`[2]
  }
  if(m==3){
    ### Test if land effect alone is different than the effect of land when pop decrease
    p_values_pop[nv]   <- linearHypothesis(model_list[[m]], test1 , vcov=vcov_list[[m]])$`Pr(>F)`[2]
    ### Test if land effect alone is different than the effect of land when le strong
    p_values_power[nv]  <- linearHypothesis(model_list[[m]], test2 , vcov=vcov_list[[m]])$`Pr(>F)`[2]
    ### Test if land effect alone is different than the effect of land when pop_decrease and le strong
    p_values_double[nv] <- linearHypothesis(model_list[[m]], test3, vcov=vcov_list[[m]])$`Pr(>F)`[2]
  }
  
### End loop1
  }
### End loop 2
}

#### Table A2
#### Tabulate 
##
pop_effect1 <- cbind.data.frame(pop_effect, p_values_pop)
pow_effect1 <- cbind.data.frame(power_effect, p_values_power)
double_effect1 <- cbind.data.frame(double_effect, p_values_double)
### 
table_wald <- cbind.data.frame(pop_effect1, pow_effect1, double_effect1)
table_wald <- round(table_wald,4)
colnames(table_wald) <- c("Pop. Effect", "P-val", "LE Power Effect", "P-val", "Pop. and Power Effect", "P-val")
table_wald[9:12,] # To show the interactions
#setwd("C://Users/franc/Dropbox/PhD/dissertation/paper_colonialism_liberalims/tex/new_version/final_tables")
#write.csv(table_wald,"table_wald.csv" )



############################################
#### Figure A2: Interaction  
#### Plot simulated coefficients from triple interaction
####
###########################################
# Dont forget to call library(Zelig)
z1 <- zelig(protect_new2 ~tierras_d*pop_decrease+ abusos_d*pop_decrease+impuestos_d*pop_decrease + 
              comunidad_d*pop_decrease+as.factor(name)+year,
            data=agn_indios_sub ,  model="logit")

z2 <- zelig(protect_new2 ~tierras_d*power+ abusos_d*power +impuestos_d*power + 
              comunidad_d*power+as.factor(name)+year,
            data=agn_indios_sub ,  model="logit")

z3 <- zelig(protect_new2 ~tierras_d*power*pop_decrease+ abusos_d*power*pop_decrease +impuestos_d*power*pop_decrease + 
              comunidad_d*power*pop_decrease+as.factor(name)+year,
            data=agn_indios_sub ,  model="logit")
#----------------------------------------------#
# Loop to create datasets 
df_pop_effect        <- data.frame()
df_le_effect         <- data.frame()
df_le_effct_pop_dec  <- data.frame()
df_le_effct_pop_inc  <- data.frame()

v <- list(c(1,0,0,0), c(0,1,0,0), c(0,0,1,0), c(0,0,0,1))
topics <- c("Land", "Mistreatment", "Taxes","Community")

#### First difference approach
### Loop across topics
for(t in 1:length(topics)){
  ### Decreasing pob
  ### Pop
  y_dec   <- setx(z1, tierras_d=v[[t]][1], abusos_d=v[[t]][2], impuestos_d=v[[t]][3], comunidad_d=v[[t]][4], pop_decrease=1)
  y_inc  <- setx(z1, tierras_d=v[[t]][1], abusos_d=v[[t]][2], impuestos_d=v[[t]][3], comunidad_d=v[[t]][4], pop_decrease=0)
  set.seed(10000)
  s.cw <- sim(z1, x = y_inc, x1 = y_dec )
  vals <- unlist(s.cw$sim.out[["x1"]][["fd"]])
  df_pop_effect <- df_pop_effect %>%  rbind.data.frame(c(mean(vals), quantile(vals, c(.025, .975))))
  
  ### Power 
  y_str   <- setx(z2, tierras_d=v[[t]][1], abusos_d=v[[t]][2], impuestos_d=v[[t]][3], comunidad_d=v[[t]][4], power = 1)
  y_weak  <- setx(z2, tierras_d=v[[t]][1], abusos_d=v[[t]][2], impuestos_d=v[[t]][3], comunidad_d=v[[t]][4], power = 0)
  set.seed(10000)
  s.cw <- sim(z2, x = y_weak, x1 = y_str  )
  vals <- unlist(s.cw$sim.out[["x1"]][["fd"]])
  df_le_effect <- df_le_effect %>% rbind.data.frame(c(mean(vals), quantile(vals, c(.025, .975))))
  
  ### Triple interaction
  #### Pop decrease and powerful LE
  y_dec_str   <- setx(z3, tierras_d=v[[t]][1], abusos_d=v[[t]][2], impuestos_d=v[[t]][3], comunidad_d=v[[t]][4], pop_decrease=1, power = 1)
  y_inc_str  <- setx(z3, tierras_d=v[[t]][1], abusos_d=v[[t]][2], impuestos_d=v[[t]][3], comunidad_d=v[[t]][4], pop_decrease=1,  power = 0)
  set.seed(10000)
  s.cw <- sim(z3, x = y_inc_str  , x1 = y_dec_str)
  vals <- unlist(s.cw$sim.out[["x1"]][["fd"]])
  df_le_effct_pop_dec <- df_le_effct_pop_dec %>% rbind.data.frame(c(mean(vals), quantile(vals, c(.025, .975))))
  
  #### Pop increase and powerful LE
  y_inc_str   <- setx(z3, tierras_d=v[[t]][1], abusos_d=v[[t]][2], impuestos_d=v[[t]][3], comunidad_d=v[[t]][4],  pop_decrease=0, power = 1)
  y_inc_weak  <- setx(z3, tierras_d=v[[t]][1], abusos_d=v[[t]][2], impuestos_d=v[[t]][3], comunidad_d=v[[t]][4],   pop_decrease=0, power = 0)
  set.seed(10000)
  s.cw <- sim(z3, x = y_inc_weak  , x1 = y_inc_str)
  vals <- unlist(s.cw$sim.out[["x1"]][["fd"]])
  df_le_effct_pop_inc <- df_le_effct_pop_inc %>% rbind.data.frame(c(mean(vals), quantile(vals, c(.025, .975))))
}

##
df_pop_effect$topic <- topics; df_le_effect$topic  <- topics
df_le_effct_pop_dec$topic <- topics ; df_le_effct_pop_inc$topic <- topics
#
df_pop_effect$effect <- "Population decrease effect" ; df_le_effect$effect  <- "Strong LE effect"
df_le_effct_pop_dec$effect <- "LE effect when pop. declining" ; df_le_effct_pop_inc$effect <- "LE effect when pop. not declining"
#
n <- c("coef", "min", "max", "topic", "Effect")
names(df_pop_effect) <- n ; names(df_le_effect) <-n ; names(df_le_effct_pop_dec) <-n ; names(df_le_effct_pop_inc)<-n 
#
df_all <- rbind.data.frame(df_pop_effect, df_le_effect, df_le_effct_pop_dec, df_le_effct_pop_inc )
df_all$order <- rep(c(1:4), 4) ; df_all$effect_order <- rep(1:4,each=4)

p <- ggplot(df_all, aes(y=coef, ymin=min, ymax=max, x=reorder(topic, -order), fill=reorder(Effect, -effect_order))) + 
         geom_point(shape=rep(c(15:18),each=4), position = position_dodge(width = 1/2), size=2.5) +  geom_errorbar(position = position_dodge(width = 1/2), width = 0) 
p <- p + geom_hline(yintercept = 0, colour = gray(1/2), lty = 2)
p <- p + theme_bw()  + ylab("First Differences") + xlab("")
p <- p  + guides(fill = guide_legend(nrow=2,byrow=TRUE, override.aes = list(shape = c(15:18)))) 
p <- p  + theme(legend.position="bottom", legend.direction='vertical', legend.title = element_blank(), legend.box.background = element_rect(colour = "black"),
                                                                                     panel.border = element_rect(colour = "black", fill=NA))
p <- p + coord_flip()
p <- p + scale_fill_discrete(labels = c("Population decline", "Strong local elites", "Strong local elites \nwith population decline", 
                                                         "Strong local elites \nwithout population decline")) 
p
ggsave(p, file="C:/Users/franc/Dropbox/PhD/dissertation/paper_colonialism_liberalims/tex/FINAL_version_2021/figures_paper/figure_a2.tiff", width=6)
#https://stackoverflow.com/questions/67640077/how-to-rotate-the-legend-labels-to-match-the-orientation-in-the-plot

######################################################################
##
##   Merits of the Case (capitulated  cases)
##   Appendix 4 (Table A4.2)
##   The assumption is that these cases are comparable in complexity
##   Only cases that were looked at more carefully
######################################################################
#### Fojas
##### quartiles using the entire sample (repetition  is not relevant here )
quart_fojas2 <- summary(agn_indios$fojas2)[5]
hist(agn_indios$fojas2, breaks=50)
#### Descriptive table (TABLE A4.2)
fojas_1      <- table(agn_indios$fojas2==1)
fojas_2_5    <- table(agn_indios$fojas2>=2 &agn_indios$fojas2<=5  )
fojas_5_10   <- table(agn_indios$fojas2>=5 &agn_indios$fojas2<10 )
fojas_10_50  <- table(agn_indios$fojas2>=10 &agn_indios$fojas2<50 )
fojas_50_132 <- table(agn_indios$fojas2>=50 &agn_indios$fojas2<132  )
fojas_132    <- table(agn_indios$fojas2>=132 )

desc_fojas <- rbind.data.frame(c("1", fojas_1[2],fojas_1[2]/nrow(agn_indios) ), 
                               c("2-5", fojas_2_5[2],fojas_2_5[2]/nrow(agn_indios) ), 
                               c("5-10", fojas_5_10[2],fojas_5_10[2]/nrow(agn_indios) ), 
                               c("10-50", fojas_10_50[2],fojas_10_50[2]/nrow(agn_indios) ), 
                               c("50-132", fojas_50_132[2],fojas_50_132[2]/nrow(agn_indios) ), 
                               c("132>", fojas_132[2],fojas_132[2]/nrow(agn_indios) )
)

colnames(desc_fojas) <- c("Pages", "No. of Files", "%")
print(xtable(desc_fojas),include.rownames=FALSE)

### add high merit to main dataset
agn_indios_sub <- agn_indios_sub %>% mutate(high_merit=ifelse(fojas2>=132,1,0))

#### subset
agn_indios_high_fojas <- agn_indios_sub %>%filter(fojas2>=quart_fojas2)

########################################
### Table A4.4: High merit as a dependent variable (to asses selection bias)
######################################

model_merit1 <- lm(high_merit~pob_ln+pop_decrease*power+ as.factor(name)*year+
                     as.factor(year), data=agn_indios_sub)

model_merit1_vcov <- cluster.vcov(model_merit1 , ~ agn_indios_sub$name)
model_merit1_coef <- coeftest(model_merit1  ,model_merit1_vcov)

### with topics
model_merit2 <- lm(high_merit~pob_ln+pop_decrease*power+tierras_d + abusos_d + impuestos_d+ comunidad_d+as.factor(name)*year+
                     as.factor(year), data=agn_indios_sub)

model_merit2_vcov <- cluster.vcov(model_merit2 , ~ agn_indios_sub$name)
model_merit2_coef <- coeftest(model_merit2  ,model_merit2_vcov)


#### Table A4.4 ####
stargazer(model_merit1 ,model_merit2, 
          se   = list(model_merit1_coef[,2], model_merit2_coef[,2]
          ) ,
          p    = list(model_merit1_coef[,4], model_merit2_coef[,4]
          ),
          dep.var.labels= c("High Complexity"),
          covariate.labels =c("Pop. Decline", "Strong LE",  "Land Conflicts", "Mistreatment", "Taxes", "Community","Pop. Decline x Strong LE"
          ),
          column.separate = c(1,1),
          object.names = F,
          #model.names = c("", "Decrease", "Increase", "Strong", "Weak", "Strong", "Weak"),
          keep = c("pop_decrease", "power", "pop_decrease:power", "tierras_d", "abusos_d", "impuestos_d", "comunidad_d"), 
          keep.stat = c("n", "adj.rsq"),
          add.lines=list(c("Province FE", "Y", "Y"),
                         c("Time FE", "Y", "Y"), 
                         c("Province time trends", "Y", "Y")
          ),
          label="high_meritt",
          title =("High Merit and Local Balance of Powers")
          #out="model_vulner_merit.tex"
)



######################################################
###  Table Table A4.5:: High Merit Subsample 
### 
######################################################
## Province and Fixed effecs #####################
#### All sample
model_dict_all_f      <- lm(protect_new2 ~tierras_d+ abusos_d +impuestos_d + comunidad_d+as.factor(year)+name*year,
                       data=agn_indios_high_fojas)
model_dict_all_f_vcov <- cluster.vcov(model_dict_all_f, ~ agn_indios_high_fojas$name)
model_dict_all_f_vcov <- coeftest(model_dict_all_f , model_dict_all_f_vcov)

#### Population variation
####  Population Decrease
model_dict_pop_dec_f <- lm(protect_new2 ~tierras_d+ abusos_d +impuestos_d + comunidad_d+ as.factor(decade)+name*decade,
                           data=agn_indios_high_fojas[which(agn_indios_high_fojas$pop_ratio<1),])
model_dict_pop_dec_f_cov <- cluster.vcov(model_dict_pop_dec_f , ~ agn_indios_high_fojas[which(agn_indios_high_fojas$pop_ratio<1),]$name)
model_dict_pop_dec_f_vcov <- coeftest(model_dict_pop_dec_f ,model_dict_pop_dec_f_cov)


#### Population Increase 
model_dict_pop_f_inc <- lm(protect_new2 ~tierras_d+ abusos_d +impuestos_d + comunidad_d + as.factor(decade)+name*decade,
                           data=agn_indios_high_fojas[which(agn_indios_high_fojas$pop_ratio>=1
                           ),])

model_dict_pop_inc_f_vcov <- cluster.vcov(model_dict_pop_f_inc , ~ agn_indios_high_fojas[which(agn_indios_high_fojas$pop_ratio>=1),]$name)
model_dict_pop_inc_f_vcov <- coeftest(model_dict_pop_f_inc , model_dict_pop_inc_f_vcov)

#### LE Power Variation
#### LE Strong
#### LE Strong
model_dict_le_f_strong <- lm(protect_new2 ~tierras_d+ abusos_d +impuestos_d +comunidad_d+as.factor(year)+name*year,
                             data=agn_indios_high_fojas[which(agn_indios_high_fojas$per_encomienda>=enc_mean
                             ),])

model_dict_le_strong_f_vcov <- cluster.vcov(model_dict_le_f_strong , ~ agn_indios_high_fojas[which(agn_indios_high_fojas$per_encomienda>=enc_mean),]$name)
model_dict_le_strong_f_vcov <- coeftest(model_dict_le_f_strong, model_dict_le_strong_f_vcov)


### LE weak
model_dict_le_f_weak <- lm(protect_new2 ~tierras_d+ abusos_d +impuestos_d +comunidad_d+as.factor(year)+name*year,
                           data=agn_indios_high_fojas[which(agn_indios_high_fojas$per_encomienda<enc_mean
                           ),])

model_dict_le_weak_f_vcov <- cluster.vcov(model_dict_le_f_weak , ~ agn_indios_high_fojas[which(agn_indios_high_fojas$per_encomienda<enc_mean),]$name)
model_dict_le_weak_f_vcov <- coeftest(model_dict_le_f_weak  ,model_dict_le_weak_f_vcov)

##### Balance of powers
#### Pop decreasing and...
#### LE Strong
model_dict_le_f_dec_s <- lm(protect_new2 ~tierras_d+ abusos_d +impuestos_d +comunidad_d+as.factor(year)+name*year,
                            data=agn_indios_high_fojas[which(agn_indios_high_fojas$per_encomienda>=enc_mean &
                                                          agn_indios_high_fojas$pop_ratio<1
                            ),])

model_dict_le_dec_s_f_vcov <- cluster.vcov(model_dict_le_f_dec_s  , ~ agn_indios_high_fojas[which(agn_indios_high_fojas$per_encomienda>=enc_mean &
                                                                                               agn_indios_high_fojas$pop_ratio<1),]$name)
model_dict_le_dec_s_f_vcov <- coeftest(model_dict_le_f_dec_s   ,model_dict_le_dec_s_f_vcov)

### LE weak
model_dict_le_f_dec_w <- lm(protect_new2 ~tierras_d+ abusos_d +impuestos_d +comunidad_d+as.factor(year)+name*year,
                            data=agn_indios_high_fojas[which(agn_indios_high_fojas$per_encomienda<enc_mean &
                                                          agn_indios_high_fojas$pop_ratio<1
                            ),])

model_dict_le_dec_w_f_vcov <- cluster.vcov(model_dict_le_f_dec_w  , ~ agn_indios_high_fojas[which(agn_indios_high_fojas$per_encomienda<enc_mean &
                                                                                               agn_indios_high_fojas$pop_ratio<1),]$name)
model_dict_le_dec_w_f_vcov <- coeftest(model_dict_le_f_dec_w   ,model_dict_le_dec_w_f_vcov)

##### Number of Units by regression
all   <- length(unique(agn_indios_high_fojas$name[which(agn_indios_high_fojas$year>1592 &
                                               agn_indios_high_fojas$name!="Mexico" )]))

dec  <- length(unique(agn_indios_high_fojas$name[which(agn_indios_high_fojas$year>1592 &
                                              agn_indios_high_fojas$name!="Mexico" & agn_indios_high_fojas$pop_ratio<1 )]))

inc  <- length(unique(agn_indios_high_fojas$name[which(agn_indios_high_fojas$year>1592 &
                                              agn_indios_high_fojas$name!="Mexico" & agn_indios_high_fojas$pop_ratio>=1 )]))

strg  <- length(unique(agn_indios_high_fojas$name[which(agn_indios_high_fojas$year>1592 &
                                               agn_indios_high_fojas$name!="Mexico" & agn_indios_high_fojas$per_encomienda>=enc_mean )]))

weak <- length(unique(agn_indios_high_fojas$name[which(agn_indios_high_fojas$year>1592 &
                                              agn_indios_high_fojas$name!="Mexico" & agn_indios_high_fojas$per_encomienda<enc_mean )]))

dec_str <- length(unique(agn_indios_high_fojas$name[which(agn_indios_high_fojas$year>1592 &
                                                 agn_indios_high_fojas$name!="Mexico" & agn_indios_high_fojas$pop_ratio<1 & agn_indios_high_fojas$per_encomienda>=enc_mean )]))

dec_weak <- length(unique(agn_indios_high_fojas$name[which(agn_indios_high_fojas$year>1592 &
                                                  agn_indios_high_fojas$name!="Mexico" & agn_indios_high_fojas$pop_ratio<1 & agn_indios_high_fojas$per_encomienda<enc_mean)]))

unit_line <- c("Units", all, dec, inc, strg, weak, dec_str, dec_weak)

#### Table A4.5 #####
#### With Robust errors
stargazer(model_dict_all_f , model_dict_pop_dec_f, model_dict_pop_f_inc, model_dict_le_f_strong, model_dict_le_f_weak,
          model_dict_le_f_dec_s, model_dict_le_f_dec_w ,
          se   = list(model_dict_all_f_vcov[, 2], model_dict_pop_dec_f_vcov[, 2], 
                      model_dict_pop_inc_f_vcov[, 2], model_dict_le_strong_f_vcov[, 2],model_dict_le_weak_f_vcov[,2],
                      model_dict_le_dec_s_f_vcov[,2], model_dict_le_dec_w_f_vcov[,2]
          ) ,
          p    = list(model_dict_all_f_vcov[, 4], model_dict_pop_dec_f_vcov[, 4], 
                      model_dict_pop_inc_f_vcov[, 4], model_dict_le_strong_f_vcov[, 4],model_dict_le_weak_f_vcov[,4],
                      model_dict_le_dec_s_f_vcov[,4], model_dict_le_dec_w_f_vcov[,4]
          ),
          dep.var.labels= c("Favorable Ruling"),
          covariate.labels =c("Land Conflicts","Mistreatment",  "Taxes", 
                              "Community"
          ),
          column.labels = c("All", "Population", "Local Elites", "Decreasing Population"
          ),
          column.separate = c(1,2,2,2),
          object.names = F,
          #model.names = c("", "Decrease", "Increase", "Strong", "Weak", "Strong", "Weak"),
          keep = c("tierras_d", "abusos_d", "impuestos_d" , "comunidad_d"), 
          keep.stat = c("n", "adj.rsq"),
          add.lines=list(c("Province FE", "Y", "Y", "Y", "Y", "Y", "Y","Y"),
                         c("Time FE", "Y", "Y", "Y", "Y","Y", "Y","Y"), 
                         c("Province time trends", "Y", "Y", "Y", "Y","Y", "Y","Y"), 
                         unit_line
          ),
          label="dictionary_vulner_merit",
          title =("Favorable Ruling, Topics, and Balance of Powers(High Merit Cases)") 
          #out="model_vulner_merit.tex"
          )



#########---------------------------------------------
######### Calculate rates of protection and merit 
######## Appendix S8
####Protection, merit and LE

## Total rate of protection
table(agn_indios_sub$protect_new2)
### Rate of protection considering power
table(agn_indios_sub$power, agn_indios_sub$protect_new2)[2,]/(sum(table(agn_indios_sub$power, agn_indios_sub$protect_new2)[2,]))
### Rate of protection considering merit
table(agn_indios_sub$high_merit, agn_indios_sub$protect_new2)[,2]/sum(table(agn_indios_sub$high_merit, agn_indios_sub$protect_new2)[,2])
### Rate of protection considering power and merit 
# High and weak
table(agn_indios_sub$power[agn_indios_sub$high_merit==1], agn_indios_sub$protect_new2[agn_indios_sub$high_merit==1])[1,]/
  sum(table(agn_indios_sub$power[agn_indios_sub$high_merit==1], agn_indios_sub$protect_new2[agn_indios_sub$high_merit==1])[1,])
# High and strong
table(agn_indios_sub$power[agn_indios_sub$high_merit==1], agn_indios_sub$protect_new2[agn_indios_sub$high_merit==1])[2,]/
sum(table(agn_indios_sub$power[agn_indios_sub$high_merit==1], agn_indios_sub$protect_new2[agn_indios_sub$high_merit==1])[2,])

### Rate of protection considering power and merit
## Low and weak
table(agn_indios_sub$power[agn_indios_sub$high_merit==0], agn_indios_sub$protect_new2[agn_indios_sub$high_merit==0])[1,]/
  sum(table(agn_indios_sub$power[agn_indios_sub$high_merit==0], agn_indios_sub$protect_new2[agn_indios_sub$high_merit==0])[1,])
# Low and strong
table(agn_indios_sub$power[agn_indios_sub$high_merit==0], agn_indios_sub$protect_new2[agn_indios_sub$high_merit==0])[2,]/
sum(table(agn_indios_sub$power[agn_indios_sub$high_merit==0], agn_indios_sub$protect_new2[agn_indios_sub$high_merit==0])[2,])

####Protection, merit and decrease
table(agn_indios_sub$pop_decrease, agn_indios_sub$high_merit)[2,]/
sum(table(agn_indios_sub$pop_decrease, agn_indios_sub$high_merit)[2,]) 
### Rate of protection considering decrease
table(agn_indios_sub$pop_decrease, agn_indios_sub$protect_new2)[2,]/
sum(table(agn_indios_sub$pop_decrease, agn_indios_sub$protect_new2)[2,])
### Rate of protection considering power and merit
# Low 
table(agn_indios_sub$pop_decrease[agn_indios_sub$high_merit==1], agn_indios_sub$protect_new2[agn_indios_sub$high_merit==1])[1,]/
  sum(table(agn_indios_sub$pop_decrease[agn_indios_sub$high_merit==1], agn_indios_sub$protect_new2[agn_indios_sub$high_merit==1])[1,])
# High
table(agn_indios_sub$pop_decrease[agn_indios_sub$high_merit==1], agn_indios_sub$protect_new2[agn_indios_sub$high_merit==1])[2,]/
sum(table(agn_indios_sub$pop_decrease[agn_indios_sub$high_merit==1], agn_indios_sub$protect_new2[agn_indios_sub$high_merit==1])[2,])

### Rate of protection considering power and merit (low)
#Low
table(agn_indios_sub$pop_decrease[agn_indios_sub$high_merit==0], agn_indios_sub$protect_new2[agn_indios_sub$high_merit==0])[1,]/
sum(table(agn_indios_sub$pop_decrease[agn_indios_sub$high_merit==0], agn_indios_sub$protect_new2[agn_indios_sub$high_merit==0])[1,])  
# High
table(agn_indios_sub$pop_decrease[agn_indios_sub$high_merit==0], agn_indios_sub$protect_new2[agn_indios_sub$high_merit==0])[2,]/
  sum(table(agn_indios_sub$pop_decrease[agn_indios_sub$high_merit==0], agn_indios_sub$protect_new2[agn_indios_sub$high_merit==0])[2,])  



######################################################################
##
##   Political Effects (Supplementary Materials S1)
##   Tables  S1.1, S1.2, S1.3
##   Figure  S1.2 and S1.3
######################################################################

#############
## Figure S1.2
## Total cases across time 
############
agn_year <- agn_indios %>% group_by(year) %>% filter(year>=1600, year<=1820) %>% summarize(total=n())

g <- ggplot(agn_year, aes(x=year, y=total)) + geom_line() + geom_smooth(span=0.3)
g <- g + theme_bw() + ylim(0,375)
g 

#ggsave(g, file="C:/Users/franc/Dropbox/PhD/dissertation/paper_colonialism_liberalims/tex/2019_version/online_appendix/total_files.png", width=8)

#################
## Test results before and after 1700 (Bourbon Reform)
## Table S1.1
#### Before 
model_dict_before1700 <- lm(protect_new2 ~tierras_d+ abusos_d +impuestos_d + comunidad_d+as.factor(year)+name*year,
                            data=agn_indios_sub[which(agn_indios_sub$year>1592 & agn_indios_sub$year<1700   ),])

model_dict_before1700_vcov <- cluster.vcov(model_dict_before1700 , ~ agn_indios_sub[which(agn_indios_sub$year>1592 & agn_indios_sub$year<1700),]$name)
model_dict_before1700_vcov <- coeftest(model_dict_before1700 , model_dict_before1700_vcov)

# After 
model_dict_after1700<- lm(protect_new2 ~tierras_d+ abusos_d +impuestos_d + comunidad_d+as.factor(year)+name*year,
                          data=agn_indios_sub[which(agn_indios_sub$year>=1700 & agn_indios_sub$year<1821  ),])
model_dict_after1700_vcov <- cluster.vcov(model_dict_after1700 , ~ agn_indios_sub[which(agn_indios_sub$year>=1700 & agn_indios_sub$year<1821 ),]$name)
model_dict_after1700_vcov <- coeftest(model_dict_after1700 , model_dict_after1700_vcov)

#################  Conditioned by period after Bourbon reform only 
# After 1700 decrease pop
model_dict_after1700_dec<- lm(protect_new2 ~tierras_d+ abusos_d +impuestos_d + comunidad_d+as.factor(decade)+name,
                              data=agn_indios_sub[which(agn_indios_sub$year>=1700 & agn_indios_sub$year<1821 & agn_indios_sub$pop_ratio<1
                              ),])
model_dict_after1700_dec_vcov <- cluster.vcov(model_dict_after1700_dec, ~ agn_indios_sub[which(agn_indios_sub$year>=1700 & agn_indios_sub$year<1821 &
                                                                                              agn_indios_sub$pop_ratio<1),]$name)
model_dict_after1700_dec_vcov <- coeftest(model_dict_after1700_dec, model_dict_after1700_dec_vcov)

# After 1700 increase pop
model_dict_after1700_inc <- lm(protect_new2 ~tierras_d+ abusos_d +impuestos_d + comunidad_d+as.factor(decade)+name,
                               data=agn_indios_sub[which(agn_indios_sub$year>=1700 & agn_indios_sub$year<1821 &
                                                        agn_indios_sub$pop_ratio>=1
                               ),])
model_dict_after1700_inc_vcov <- cluster.vcov(model_dict_after1700_inc, ~ agn_indios_sub[which(agn_indios_sub$year>=1700 & agn_indios_sub$year<1821 &
                                                                                              agn_indios_sub$pop_ratio>=1),]$name)
model_dict_after1700_inc_vcov <- coeftest(model_dict_after1700_inc, model_dict_after1700_inc_vcov)

# After 1700 le strong
## Recalculate encomienda mean for this period
enc_mean_1700 <- mean(agn_indios_sub$per_encomienda[agn_indios_sub$year>=1700 & agn_indios_sub$year<1821],na.rm=T)

model_dict_after1700_str <- lm(protect_new2 ~tierras_d+ abusos_d +impuestos_d + comunidad_d+as.factor(decade)+name,
                               data=agn_indios_sub[which(agn_indios_sub$year>=1700 & agn_indios_sub$year<1821 &
                                                       agn_indios_sub$per_encomienda>=enc_mean_1700
                               ),])
model_dict_after1700_str_vcov <- cluster.vcov(model_dict_after1700_str, ~ agn_indios_sub[which(agn_indios_sub$year>=1700 & agn_indios_sub$year<1821 &
                                                                                             agn_indios_sub$per_encomienda>=enc_mean_1700),]$name)
model_dict_after1700_str_vcov <- coeftest(model_dict_after1700_str, model_dict_after1700_str_vcov)

# After 1700 le weak
model_dict_after1700_wk <- lm(protect_new2 ~tierras_d+ abusos_d +impuestos_d + comunidad_d+as.factor(decade)+name,
                              data=agn_indios_sub[which(agn_indios_sub$year>=1700 & agn_indios_sub$year<1821 &
                                                       agn_indios_sub$per_encomienda<enc_mean_1700
                              ),])
model_dict_after1700_wk_vcov <- cluster.vcov(model_dict_after1700_wk, ~ agn_indios_sub[which(agn_indios_sub$year>=1700 & agn_indios_sub$year<1821 &
                                                                                           agn_indios_sub$per_encomienda<enc_mean_1700),]$name)
model_dict_after1700_wk_vcov <- coeftest(model_dict_after1700_wk, model_dict_after1700_wk_vcov )

##### Calculate # of Units 
bef <- length(unique(agn_indios_sub$name[agn_indios_sub$year>1592 & agn_indios_sub$year<1700]))
aft <- length(unique(agn_indios_sub$name[agn_indios_sub$year>=1700 & agn_indios_sub$year<1821 ]))
## Conditioned on after
aft_d <- length(unique(agn_indios_sub$name[agn_indios_sub$year>=1700 & agn_indios_sub$year<1821  & agn_indios_sub$pop_ratio<1]))
aft_i <- length(unique(agn_indios_sub$name[agn_indios_sub$year>=1700 & agn_indios_sub$year<1821 & agn_indios_sub$pop_ratio>=1]))
aft_s <- length(unique(agn_indios_sub$name[agn_indios_sub$year>=1700 & agn_indios_sub$year<1821 & agn_indios_sub$per_encomienda<enc_mean_1700]))
aft_w <- length(unique(agn_indios_sub$name[agn_indios_sub$year>=1700 & agn_indios_sub$year<1821 & agn_indios_sub$per_encomienda>=enc_mean_1700]))
aft_d_s <- length(unique(agn_indios_sub$name[agn_indios_sub$year>=1700 & agn_indios_sub$year<1821 & agn_indios_sub$pop_ratio<1 &
                                           agn_indios_sub$per_encomienda<enc_mean_1700]))
aft_d_w <- length(unique(agn_indios_sub$name[agn_indios_sub$year>=1700 & agn_indios_sub$year<1821 & agn_indios_sub$pop_ratio>=1
                                         & agn_indios_sub$per_encomienda>=enc_mean_1700
                                         ]))


### Table S1.1 (Supplementary materials)
stargazer(model_dict_before1700, model_dict_after1700, model_dict_after1700_dec, model_dict_after1700_inc ,model_dict_after1700_str, model_dict_after1700_wk,
          se   = list(model_dict_before1700_vcov[,2], model_dict_after1700_vcov[,2],
                      model_dict_after1700_dec_vcov[,2], model_dict_after1700_inc_vcov[,2] ,model_dict_after1700_str_vcov[,2], 
                      model_dict_after1700_wk_vcov[,2]
          ) ,
          p    = list(model_dict_before1700_vcov[,4], model_dict_after1700_vcov[,4],
                      model_dict_after1700_dec_vcov[,4], model_dict_after1700_inc_vcov[,4] ,model_dict_after1700_str_vcov[,4], 
                      model_dict_after1700_wk_vcov[,4]
          ),
          
          dep.var.labels= c("Favorable Ruling"),
          covariate.labels =c("Land Conflicts","Mistreatment",  "Taxes", 
                              "Community"
          ),
          column.labels = c("Pre-1700", "Post-1700", "Population", "Local Elites"
          ),
          column.separate = c(1,1,2,2),
          #model.names = c("", "Decrease", "Increase", "Strong", "Weak", "Strong", "Weak"),
          keep = c("tierras_d", "abusos_d", "impuestos_d" , "comunidad_d"), 
          keep.stat = c("n", "adj.rsq", "ll"),
          add.lines=list(c("Province FE", "Y", "Y", "Y", "Y", "Y"),
                         c("Time FE", "Y", "Y", "Y", "Y", "Y"),
                         c("Province linear trend", "Y", "Y", "Y", "Y", "Y"),
                         c("Units", aft_d, aft_i, aft_s, aft_w, aft_d_s, aft_d_w)
          ),
          label="after_1700",
          title =("Indigenous Vulnerability after Bourbon reform") 
          #out="1700_after.tex"
          )


#############################################
##### Figure S1.3:Convergence effects 
#####  Coeff plot at different moments
## pre 1700
coefs_l <- model_dict_before1700_vcov[2,(1:2)]
coefs_a <- model_dict_before1700_vcov[3,(1:2)]
coefs_i <- model_dict_before1700_vcov[4,(1:2)]
coefs_c <- model_dict_before1700_vcov[5,(1:2)]
## 1700-1725
model_dict_p1 <- lm(protect_new2 ~tierras_d+ abusos_d +impuestos_d + comunidad_d+as.factor(year)+name*year,
                          data=agn_indios_sub[which(agn_indios_sub$year>=1700 & agn_indios_sub$year<1725 &
                                                      agn_indios_sub$name!="Mexico"
                          ),])
model_dict_p1_vcov <- cluster.vcov(model_dict_p1 , ~ agn_indios_sub[which(agn_indios_sub$year>=1700 & agn_indios_sub$year<1725 &
                                                                                          agn_indios_sub$name!="Mexico"),]$name)
model_dict_p1_vcov <- coeftest(model_dict_p1 , model_dict_p1_vcov)

#
coefs_l  <- rbind.data.frame(coefs_l, model_dict_p1_vcov[2,(1:2)])
coefs_a  <- rbind.data.frame(coefs_a, model_dict_p1_vcov[3,(1:2)])
coefs_i  <- rbind.data.frame(coefs_i, model_dict_p1_vcov[4,(1:2)])
coefs_c  <- rbind.data.frame(coefs_c, model_dict_p1_vcov[5,(1:2)])

## 1725-1750
model_dict_p2 <- lm(protect_new2 ~tierras_d+ abusos_d +impuestos_d + comunidad_d+as.factor(year)+name*year,
                    data=agn_indios_sub[which(agn_indios_sub$year>=1725 & agn_indios_sub$year<1750 &
                                                agn_indios_sub$name!="Mexico"
                    ),])
model_dict_p2_vcov <- cluster.vcov(model_dict_p2 , ~ agn_indios_sub[which(agn_indios_sub$year>=1725 & agn_indios_sub$year<1750 &
                                                                            agn_indios_sub$name!="Mexico"),]$name)
model_dict_p2_vcov <- coeftest(model_dict_p2 , model_dict_p2_vcov)

#
coefs_l  <- rbind.data.frame(coefs_l, model_dict_p2_vcov[2,(1:2)])
coefs_a  <- rbind.data.frame(coefs_a, model_dict_p2_vcov[3,(1:2)])
coefs_i  <- rbind.data.frame(coefs_i, model_dict_p2_vcov[4,(1:2)])
coefs_c  <- rbind.data.frame(coefs_c, model_dict_p2_vcov[5,(1:2)])
## 1750-1821
model_dict_p3 <- lm(protect_new2 ~tierras_d+ abusos_d +impuestos_d + comunidad_d+as.factor(year)+name*year,
                    data=agn_indios_sub[which(agn_indios_sub$year>=1750 & agn_indios_sub$year<1821 &
                                                agn_indios_sub$name!="Mexico"
                    ),])
model_dict_p3_vcov <- cluster.vcov(model_dict_p3 , ~ agn_indios_sub[which(agn_indios_sub$year>=1750 & agn_indios_sub$year<1821 &
                                                                            agn_indios_sub$name!="Mexico"),]$name)
model_dict_p3_vcov <- coeftest(model_dict_p3 , model_dict_p3_vcov)

#
coefs_l  <- rbind.data.frame(coefs_l, model_dict_p3_vcov[2,(1:2)])
coefs_a  <- rbind.data.frame(coefs_a, model_dict_p3_vcov[3,(1:2)]) 
coefs_i  <- rbind.data.frame(coefs_i, model_dict_p3_vcov[4,(1:2)])
coefs_c  <- rbind.data.frame(coefs_c, model_dict_p3_vcov[5,(1:2)])
####
colnames(coefs_l) <- c("coef", "se") ; colnames(coefs_a) <- c("coef", "se");colnames(coefs_i) <- c("coef", "se");colnames(coefs_c) <- c("coef", "se")
periods <- c("Pre-1700", "1700-1725", "1725-1750","1750-1821")
coefs_l$periods <- periods ; coefs_a$periods <- periods; coefs_i$periods <- periods ; coefs_c$periods <- periods
coefs_l$order <- c(1:4) ; coefs_a$order <- c(1:4) ; coefs_i$order <- c(1:4)  ; coefs_c$order <- c(1:4)
### ggplots
## land
l <- ggplot(coefs_l, aes(x=reorder(periods, order), y=coef, ymin=coef-1.96*se, ymax=coef+1.96*se, group=1)) + geom_pointrange() + geom_line()
l <- l + theme_bw() + ylab("Coefficient") + xlab("") + geom_hline(yintercept = 0, color="red", linetype="dashed")
l <- l + ggtitle("Land") + theme(axis.text.x = element_text(angle = 90)) 
l
## abuses
a <- ggplot(coefs_a, aes(x=reorder(periods, order), y=coef, ymin=coef-1.96*se, ymax=coef+1.96*se, group=1)) + geom_pointrange() + geom_line()
a <- a + theme_bw() + ylab("") + xlab("") + geom_hline(yintercept = 0, color="red", linetype="dashed")
a <- a + ggtitle("Mistreatment")+ theme(axis.text.x = element_text(angle = 90)) 
a
## taxes
i <- ggplot(coefs_i, aes(x=reorder(periods, order), y=coef, ymin=coef-1.96*se, ymax=coef+1.96*se, group=1)) + geom_pointrange() + geom_line()
i <- i + theme_bw() + ylab("Coefficient") + xlab("") + geom_hline(yintercept = 0, color="red", linetype="dashed")
i <- i  + ggtitle("Taxes")+ theme(axis.text.x = element_text(angle = 90)) 
i
## community
c <- ggplot(coefs_c, aes(x=reorder(periods, order), y=coef, ymin=coef-1.96*se, ymax=coef+1.96*se, group=1)) + geom_pointrange() + geom_line()
c <- c + theme_bw() + ylab("") + xlab("") + geom_hline(yintercept = 0, color="red", linetype="dashed")
c <- c + ggtitle("Community")+ theme(axis.text.x = element_text(angle = 90)) 
c
#### Save grid 
# ggsave(file= "C:/Users/franc/Dropbox/PhD/dissertation/paper_colonialism_liberalims/tex/2019_version/online_appendix/bourbon.png",
#        grid.arrange(l,a,i,c, nrow = 2), width=7)


###################################
## Table S1.2
## Testing for additional political factors
######################################

##   Read dataset containing data on viceroys, kings, and foreign war involvement
viceroy <- read.csv("C:/Users/franc/Dropbox/PhD/Edgar_Tlaxcala/kings_viceroys/kings_viceroys_new_spain.csv") %>% dplyr::rename(year=Year)
viceroy <- viceroy %>% mutate_if( is.numeric, ~replace(., is.na(.), 0))

## Create a variable measuring total Spanish involvement in wars in Europe 
war <- rowSums(viceroy[5:65])
viceroy$wars <- war
## War dummy
viceroy <- viceroy %>% mutate(war_dum=ifelse(wars>0,1,0), war_dum2=ifelse(wars>=2,1,0))

#### Merge with main dataset
agn_indios_sub <- agn_indios_sub %>% left_join(viceroy )

#### All sample (unconditional)
model_dict_all <- lm(protect_new2 ~tierras_d+ abusos_d +impuestos_d + comunidad_d+as.factor(year)+name*year,
                     data=agn_indios_sub)
model_dict_all_vcov <- cluster.vcov(model_dict_all   , ~ agn_indios_sub$name)
model_dict_all_vcov <- coeftest(model_dict_all , model_dict_all_vcov)

##### All Political variables 
model_dict_politics <- lm(protect_new2 ~tierras_d+ abusos_d +impuestos_d + comunidad_d+Viceroy+King+ war_dum +as.factor(year)+name*year,
                          data=agn_indios_sub)
model_dict_politics_vcov <- cluster.vcov(model_dict_politics  , ~ agn_indios_sub$name)
model_dict_politics_vcov <- coeftest(model_dict_politics  , model_dict_politics_vcov)

####### Table A10
stargazer( model_dict_all, model_dict_politics,
           se   = list(model_dict_all_vcov[,2], model_dict_politics_vcov[,2]
           ) ,
           p    = list(model_dict_all_vcov[,4],  model_dict_politics_vcov[,4]
           ),
           dep.var.labels= c("Favorable Ruling"),
           covariate.labels =c("Land Conflicts","Mistreatment",  "Taxes", 
                               "Community"
           ),
           column.labels = c("Baseline Model", "Internal and External Politics"
           ),
           column.separate = c(1,1),
           object.names = F,
           #model.names = c("", "Decrease", "Increase", "Strong", "Weak", "Strong", "Weak"),
           keep = c("tierras_d", "abusos_d", "impuestos_d" , "comunidad_d"), 
           keep.stat = c("n", "adj.rsq"),
           add.lines=list(c("Province controls", "Y", "Y"),
                          c("Time controls", "Y", "Y"),
                          c("Province time trend", "Y", "Y"),
                          c("Viceroy FE", "N", "Y"),
                          c("King FE", "N", "Y"),
                          c("War FE", "N", "Y")
           ),
           label="dictionary_vulner_merit",
           title =("Favorable Ruling and Political Factors)")
           #out="model_politics.tex"
           )

####################################
### Table S1.3
### Testing direct vs Indirect Ruling  
####################################

## This files contains population and percentage under encomienda data by decade (not interpolated)
panel_pob_decades <- read.csv("C://Users/franc/Dropbox/PhD/Edgar_Tlaxcala/AGN/data/pob_decades.csv") #In Tlaxcala agn 
panel_pob_decades <- panel_pob_decades %>% dplyr::select(name, decade, decade_id, pob, per_encomienda) 

# Subset of provinces with no encomienda (places that never had encomienda or dissaperaed almost immediately)
no_enc     <- panel_pob_decades %>% filter(decade==1580 & per_encomienda==0) %>% mutate(ever_encomienda=0) %>% dplyr::select(name, ever_encomienda)

#### Merge and create ever_encomienda
agn_indios_sub <- agn_indios_sub %>% left_join(no_enc) %>% mutate(ever_enc=ifelse(ever_encomienda==0 & !is.na(per_encomienda),0,1))

#### Provinces than never had an encomienda
model_dict_no_enc <- lm(protect_new2 ~tierras_d+ abusos_d +impuestos_d + comunidad_d+as.factor(year)+name*year,
                        data=agn_indios_sub[which(agn_indios_sub$year>1592 & agn_indios_sub$ever_enc==0
                        ),])
model_dict_no_enc_vcov  <- cluster.vcov(model_dict_no_enc   , ~ agn_indios_sub[which(agn_indios_sub$year>1592  &  agn_indios_sub$ever_enc==0),]$name)
model_dict_no_enc_vcov  <- coeftest(model_dict_no_enc , model_dict_no_enc_vcov )


#### Provinces that had encomienda at any point in time 
model_dict_enc <- lm(protect_new2 ~tierras_d+ abusos_d +impuestos_d + comunidad_d+as.factor(year)+name*year,
                     data=agn_indios_sub[which(agn_indios_sub$year>1592 & agn_indios_sub$ever_enc==1
                     ),])
model_dict_enc_vcov  <- cluster.vcov(model_dict_enc   , ~ agn_indios_sub[which(agn_indios_sub$year>1592 &  agn_indios_sub$ever_enc==1),]$name)
model_dict_enc_vcov  <- coeftest(model_dict_enc , model_dict_enc_vcov )

### Units
no_enc_units <- length(unique(agn_indios_sub$name[agn_indios_sub$year>1592 &  agn_indios_sub$ever_enc==0]))

yes_enc_units <- length(unique(agn_indios_sub$name[agn_indios_sub$year>1592 &   agn_indios_sub$ever_enc==1]))

##### Table A11. Direct and Indirect Ruling 
stargazer(model_dict_no_enc , model_dict_enc , 
          se   = list(model_dict_no_enc_vcov[,2] , model_dict_enc_vcov[,2]  
          ) ,
          p    = list(model_dict_no_enc_vcov[,4] , model_dict_enc_vcov[,4]
          ),
          dep.var.labels= c("Favorable Ruling"),
          covariate.labels =c("Land Conflicts","Mistreatment",  "Taxes", 
                              "Community"
          ),
          column.labels = c("Always Direct Ruling", "Indirect Ruling"),
          object.names = F,
          #model.names = c("", "Decrease", "Increase", "Strong", "Weak", "Strong", "Weak"),
          keep = c("tierras_d", "abusos_d", "impuestos_d" , "comunidad_d"), 
          keep.stat = c("n", "adj.rsq"),
          add.lines=list(c("Province controls", "Y", "Y"),
                         c("Time controls", "Y", "Y"), 
                         c("Province linear trend", "Y", "Y"),
                         c("Units", no_enc_units , yes_enc_units )
          ),
          label="dictionary_rule",
          title =("Favorable Ruling, Topics, and Balance of Powers (Direct vs. Indirect Ruling)") 
#          out="direct_indirect.tex"
)



#---------------------------------------------------#
##################################
#### Plaintiff selection bias (Supplemantary materials S9)
##################################

############################
#### Figure S9.1 
### How many cases  are favorable by topic (percentage)
##########################

land <- agn_indios %>% filter(year>=1600, year<=1820) %>% dplyr::select(year, decade, name, tierras_d, protect_new2) %>% filter(tierras_d==1) %>%
          group_by(decade) %>% summarise(tot=sum(tierras_d), per=sum(protect_new2)/sum(tierras_d), Topic="Land")

tax <- agn_indios %>% filter(year>=1600, year<=1820)  %>% dplyr::select(year, decade, name,impuestos_d, protect_new2) %>% filter(impuestos_d==1) %>%
  group_by(decade) %>% summarise(per=sum(protect_new2)/sum(impuestos_d),tot=sum(impuestos_d), Topic="Taxes")

abuse <- agn_indios %>% filter(year>=1600, year<=1820) %>% dplyr::select(year, decade, name,abusos_d, protect_new2) %>% filter(abusos_d==1) %>%
  group_by(decade) %>% summarise(per=sum(protect_new2)/sum(abusos_d), tot=sum(abusos_d), Topic="Mistreatment")

community <- agn_indios %>% filter(year>=1600, year<=1820) %>% dplyr::select(year, decade, name, comunidad_d, protect_new2) %>% filter(comunidad_d==1) %>%
  group_by(decade) %>% summarise(per=sum(protect_new2)/sum(comunidad_d), tot=sum(comunidad_d), Topic="Community")

###
df <- rbind.data.frame(land, tax, abuse, community)

g <- ggplot(df, aes(x=decade, y=tot, color=Topic)) + stat_smooth(size=1, se=F, span=0.5) 
g <- g + theme_bw() + ylim(0,250)
g <- g + xlab("") + ylab("Total Cases")
g
#ggsave("C://Users/franc/Dropbox/PhD/dissertation/paper_colonialism_liberalims/tex/2019_version/plots/decade_topic.png", g, width=8)


###########################
# Topic as dependent variable with prior rate of success 
# Table S9.1
###########################
### Create rate of protection of the previous 50 years by topic
agn_roll <- agn_indios_sub %>% mutate(protect_land   = ifelse(protect_new2==1 &tierras_d==1,1,0),
                                      protect_abusos = ifelse(protect_new2==1 &abusos_d==1,1,0),
                                      protect_taxes  = ifelse(protect_new2==1 &impuestos_d==1,1,0),
                                      protect_community   = ifelse(protect_new2==1 &comunidad_d==1,1,0)
                                      )  %>% 
                                     group_by(decade, name ) %>%
                                summarise(land_tot=sum(tierras_d), protect_land=sum(protect_land),
                                          abusos_tot=sum(abusos_d), protect_abusos=sum(protect_abusos),
                                          taxes_tot=sum(impuestos_d), protect_taxes=sum(protect_taxes),
                                          community_tot=sum(comunidad_d), protect_community=sum(protect_community)
                                          )%>% 
                                     arrange(name, decade)%>% filter(!is.na(decade) & !is.na(name))
#### Rolling sums
agn_roll$land_tot_roll<- c(NA);agn_roll$abusos_tot_roll<- c(NA);agn_roll$taxes_tot_roll<- c(NA);agn_roll$community_tot_roll <- c(NA)
agn_roll$land_prot_roll<- c(NA);agn_roll$abusos_prot_roll<- c(NA);agn_roll$taxes_prot_roll<- c(NA);agn_roll$community_prot_roll <- c(NA)

##### Loop 
for(i in 5:nrow(agn_roll)){
  if(agn_roll$name[i]==agn_roll$name[i-4]){
    ### Land 
    agn_roll$land_tot_roll[i] <- sum(agn_roll$land_tot[i], agn_roll$land_tot[i-1],agn_roll$land_tot[i-2],
                                          agn_roll$land_tot[i-3], agn_roll$land_tot[i-4])
    agn_roll$land_prot_roll[i] <- sum(agn_roll$protect_land[i], agn_roll$protect_land[i-1],agn_roll$protect_land[i-2],
                                          agn_roll$protect_land[i-3], agn_roll$protect_land[i-4])
    ##Abuses
    agn_roll$abusos_tot_roll[i] <- sum(agn_roll$abusos_tot[i], agn_roll$abusos_tot[i-1],agn_roll$abusos_tot[i-2],
                                     agn_roll$abusos_tot[i-3], agn_roll$abusos_tot[i-4])
    agn_roll$abusos_prot_roll[i] <- sum(agn_roll$protect_abusos[i], agn_roll$protect_abusos[i-1],agn_roll$protect_abusos[i-2],
                                      agn_roll$protect_abusos[i-3], agn_roll$protect_abusos[i-4])
    ##Taxes
    agn_roll$taxes_tot_roll[i] <- sum(agn_roll$taxes_tot[i], agn_roll$taxes_tot[i-1],agn_roll$taxes_tot[i-2],
                                       agn_roll$taxes_tot[i-3], agn_roll$taxes_tot[i-4])
    agn_roll$taxes_prot_roll[i] <- sum(agn_roll$protect_taxes[i], agn_roll$protect_taxes[i-1],agn_roll$protect_taxes[i-2],
                                        agn_roll$protect_taxes[i-3], agn_roll$protect_taxes[i-4])
    ##Community 
    agn_roll$community_tot_roll[i] <- sum(agn_roll$community_tot[i], agn_roll$community_tot[i-1],agn_roll$community_tot[i-2],
                                      agn_roll$community_tot[i-3], agn_roll$community_tot[i-4])
    agn_roll$community_prot_roll[i] <- sum(agn_roll$protect_community[i], agn_roll$protect_community[i-1],agn_roll$protect_community[i-2],
                                       agn_roll$protect_community[i-3], agn_roll$protect_community[i-4])
  }
}

agn_roll <- agn_roll %>% mutate(per_land_prot   =(land_prot_roll/land_tot_roll ),per_land_prot=ifelse(is.nan(per_land_prot),0,per_land_prot),
                                per_abusos_prot =(abusos_prot_roll/abusos_tot_roll ),per_abusos_prot=ifelse(is.nan(per_abusos_prot),0,per_abusos_prot),
                                per_taxes_prot  =(taxes_prot_roll/taxes_tot_roll ),per_taxes_prot=ifelse(is.nan(per_taxes_prot),0,per_taxes_prot),
                                per_community_prot =(community_prot_roll/community_tot_roll ),per_community_prot=ifelse(is.nan(per_community_prot),0,per_community_prot)
                                ) %>% 
                 dplyr::select(decade, name,per_land_prot, per_abusos_prot, per_taxes_prot, per_community_prot )

agn_indios_sub <- agn_indios_sub %>% left_join(agn_roll, by=c("name", "decade"))

#########
# Test prevalence of cases by topic conditioned on certain variables
##########
### Land 
land <- lm(tierras_d ~ power + pop_decrease + per_land_prot + as.factor(decade)+name*decade,
                          data=agn_indios_sub)
land_vcov <- cluster.vcov(land , ~ agn_indios_sub$name)
land_vcov <- coeftest(land, land_vcov)

### Abusos
abusos <- lm(abusos_d ~ power + pop_decrease + per_abusos_prot + as.factor(decade)+name*decade,
           data=agn_indios_sub)
abusos_vcov <- cluster.vcov(abusos , ~ agn_indios_sub$name)
abusos_vcov <- coeftest(abusos, abusos_vcov)

### Taxes
taxes <- lm(impuestos_d ~ power + pop_decrease + per_taxes_prot + as.factor(decade)+name*decade,
             data=agn_indios_sub)
taxes_vcov <- cluster.vcov(taxes , ~ agn_indios_sub$name)
taxes_vcov <- coeftest(taxes , taxes_vcov)

### Community 
community <- lm(comunidad_d ~ power + pop_decrease + per_community_prot + as.factor(decade)+name*decade,
             data=agn_indios_sub)
community_vcov <- cluster.vcov(community  , ~ agn_indios_sub$name)
community_vcov <- coeftest(community , community_vcov)

###### Clusters
land_clus   <- length(unique(agn_indios_sub$name[!is.na(agn_indios_sub$per_land_prot)]))
abusos_clus <- length(unique(agn_indios_sub$name[!is.na(agn_indios_sub$per_abusos_prot)]))
taxes_clus <- length(unique(agn_indios_sub$name[!is.na(agn_indios_sub$per_taxes_prot)]))
community_clus <- length(unique(agn_indios_sub$name[!is.na(agn_indios_sub$per_community_prot)]))

##### Table A11. Direct and Indirect Ruling 
stargazer(land, abusos, taxes, community, 
          se   = list(land_vcov[,2] , abusos_vcov[,2] , taxes_vcov[,2], community_vcov[,2]
          ) ,
          p    = list(land_vcov[,4] , abusos_vcov[,4] , taxes_vcov[,4], community_vcov[,4]
          ),
          dep.var.labels= c("Land", "Mistreatment", "Taxes", "Community"),
          covariate.labels =c("LE Power", "Pop. decrease", "Rate of protection land", 
                              "Rate of protection abuses","Rate of protection taxes","Rate of protection community"
                              
                              
          ),
          #column.labels = c("Always Direct Ruling", "Indirect Ruling"),
          object.names = F,
          #model.names = c("", "Decrease", "Increase", "Strong", "Weak", "Strong", "Weak"),
          keep = c("power", "pop_decrease", "per_land_prot", "per_abusos_prot", 
                            "per_taxes_prot","per_community_prot" ), 
          keep.stat = c("n", "adj.rsq"),
          add.lines=list(c("Province controls", "Y", "Y",  "Y", "Y"),
                         c("Time controls", "Y", "Y" ,"Y", "Y"), 
                         c("Province linear trend", "Y", "Y",  "Y", "Y"),
                         c("Units", land_clus  , abusos_clus , taxes_clus , community_clus )
          ),
          label="topics",
          title =("Elements for Claim Topics") 
          #          out="direct_indirect.tex"
)



############################################
### Figure S9.2
### Do provinces that use the court more learn the ways of the court?
####### High inensity provinces 
#agn_panel_pob_decades

agn_provinces_high <- agn_indios %>% filter(year>=1600, year<=1820) %>% group_by(name)  %>% 
                      summarise(files=n()) %>% filter(files>=100, !is.na(name), name!="Mexico")
length(unique(agn_provinces_high$name))
agn_provinces_high <- agn_indios[agn_indios$name %in% agn_provinces_high$name,]

agn_provinces_high <- agn_provinces_high %>% filter(year>=1600, year<=1820) %>% group_by(decade) %>% 
  summarise(files=n(), protect=sum(protect_new2)) %>% mutate(protect_per=protect/files, Intensity="High")


####### low inensity provinces 
#agn_panel_pob_decades
agn_provinces_low <- agn_indios %>% filter(year>=1600, year<=1820) %>% group_by(name)  %>% 
  summarise(files=n()) %>% filter(files<100, !is.na(name), name!="Mexico")
length(unique(agn_provinces_low$name))
agn_provinces_low <- agn_indios[agn_indios$name %in% agn_provinces_low$name,]

agn_provinces_low <- agn_provinces_low %>% filter(year>=1600, year<=1820) %>% group_by(decade) %>% 
  summarise(files=n(), protect=sum(protect_new2)) %>% mutate(protect_per=protect/files, Intensity="Low")

##########
agn_provinces_intensity <- rbind.data.frame(agn_provinces_high, agn_provinces_low)

p <- ggplot(agn_provinces_intensity , aes(x=decade, y=protect_per, color=Intensity)) + geom_line()
p <- p + theme_bw() + xlab("") + ylab("Proportion of cases with favorable ruling")
p
#ggsave("C://Users/franc/Dropbox/PhD/dissertation/paper_colonialism_liberalims/tex/2019_version/plots/decade_intensity.png", p, width=8)

##### Create dataset with high merit cases
#panel_files <- agn_indios %>% group_by(name, decade) %>% filter(fojas2>=132) %>% summarise(high_merit=n())
#write.csv(panel_files, "high_merit_panel.csv" )


#############################################
##  Supplementary material S10
##  Modeling with one topic only
#####################################

#### Table S10.1
#### Count coocurrence:
table(agn_indios$only_one)
# Land
table(agn_indios$tierras_d)  # Total
table(agn_indios$tierras_d3) # Only land
table(agn_indios$tierras_d &  agn_indios$abusos_d) # with abusos
table(agn_indios$tierras_d &  agn_indios$impuestos_d)# with taxes
table(agn_indios$tierras_d &  agn_indios$comunidad_d)# with community
table(agn_indios$tierras_d &  agn_indios$abusos_d &  agn_indios$impuestos_d)# with abusos & impuestos
table(agn_indios$tierras_d &  agn_indios$abusos_d &  agn_indios$community_d)# with abusos & comunidad
## Abuses
table(agn_indios$abusos_d)  # Total
table(agn_indios$abusos_d3) # Only abusos
table(agn_indios$abusos_d &  agn_indios$impuestos_d)# with taxes
table(agn_indios$abusos_d &  agn_indios$comunidad_d)# with community
## Taxes
table(agn_indios$impuestos_d)  # Total
table(agn_indios$impuestos_d3) # Only taxes
table(agn_indios$impuestos_d &  agn_indios$comunidad_d)# with community
## Comunidad
table(agn_indios$comunidad_d)  # Total
table(agn_indios$comunidad_d3) # Only community


############# Table S10.2
## Province and Fixed effecs #####################
#### All sample (unconditional)
model_dict2_all <- lm(protect_new2 ~theme+as.factor(year)+name*year,
                     data=agn_indios_sub)
model_dict2_all_vcov <- cluster.vcov(model_dict2_all   , ~ agn_indios_sub$name)
model_dict2_all_vcov <- coeftest(model_dict2_all , model_dict2_all_vcov)

#####
#### Conditioned on Population Decrease
model_dict2_pop_dec <- lm(protect_new2 ~theme +as.factor(year)+name*year,
                         data=agn_indios_sub[which(agn_indios_sub$pop_ratio<1),])
model_dict2_pop_dec_vcov <- cluster.vcov(model_dict2_pop_dec  , ~ agn_indios_sub[which(agn_indios_sub$pop_ratio<1),]$name)
model_dict2_pop_dec_vcov <- coeftest(model_dict2_pop_dec , model_dict2_pop_dec_vcov)

#### Population Increase 
model_dict2_pop_inc <- lm(protect_new2 ~theme+ as.factor(year)+name*year,
                         data=agn_indios_sub[which(agn_indios_sub$pop_ratio>=1),])

model_dict2_pop_inc_vcov <- cluster.vcov(model_dict2_pop_inc  , ~ agn_indios_sub[which(agn_indios_sub$pop_ratio>=1),]$name)
model_dict2_pop_inc_vcov <- coeftest(model_dict2_pop_inc , model_dict2_pop_inc_vcov)


#### LE Power
#### LE Strong
## Poulation in Encomienda mean 
enc_mean <- mean(agn_indios_sub$per_encomienda, na.rm=T) # 0.1529982

model_dict2_le_strong <- lm(protect_new2 ~theme+as.factor(year)+name*year,
                           data=agn_indios_sub[which(agn_indios_sub$per_encomienda>=enc_mean),])

model_dict2_le_strong_vcov <- cluster.vcov(model_dict2_le_strong , ~ agn_indios_sub[which(agn_indios_sub$per_encomienda>=enc_mean),]$name)
model_dict2_le_strong_vcov <- coeftest(model_dict2_le_strong , model_dict2_le_strong_vcov)


### LE weak
model_dict2_le_weak <- lm(protect_new2 ~theme+as.factor(year)+name*year,
                         data=agn_indios_sub[which(agn_indios_sub$per_encomienda<enc_mean 
                         ),])
model_dict2_le_weak_vcov <- cluster.vcov(model_dict2_le_weak , ~ agn_indios_sub[which(agn_indios_sub$per_encomienda<enc_mean),]$name)
model_dict2_le_weak_vcov <- coeftest(model_dict2_le_weak  , model_dict2_le_weak_vcov)


#### Conditioned by pop Decrease AND...
## LE strong
model_dict2_pop_dec_strong <- lm(protect_new2 ~theme+ as.factor(year)+name*year,
                                data=agn_indios_sub[which(agn_indios_sub$pop_ratio<1 & agn_indios_sub$per_encomienda>=enc_mean
                                ),])

model_dict2_pop_dec_strong_vcov <- cluster.vcov(model_dict2_pop_dec_strong , ~ agn_indios_sub[which(agn_indios_sub$pop_ratio<1 & agn_indios_sub$per_encomienda>=enc_mean),]$name)
model_dict2_pop_dec_strong_vcov <- coeftest(model_dict2_pop_dec_strong , model_dict2_pop_dec_strong_vcov)


## LE weak
model_dict2_pop_dec_weak <- lm(protect_new2 ~theme+as.factor(year)+name*year,
                              data=agn_indios_sub[which( agn_indios_sub$pop_ratio<1 & agn_indios_sub$per_encomienda<enc_mean
                              ),])

model_dict2_pop_dec_weak_vcov  <- cluster.vcov(model_dict2_pop_dec_weak  , ~ agn_indios_sub[which(agn_indios_sub$pop_ratio<1 & agn_indios_sub$per_encomienda<enc_mean),]$name)
model_dict2_pop_dec_weak_vcov<- coeftest(model_dict2_pop_dec_weak  ,model_dict2_pop_dec_weak_vcov)


#
##### Number of Units by regression
all   <- length(unique(agn_indios_sub$name[which(agn_indios_sub$only_one==1)]))

dec  <- length(unique(agn_indios_sub$name[which(agn_indios_sub$pop_ratio<1 & agn_indios_sub$only_one==1)]))

inc  <- length(unique(agn_indios_sub$name[which(agn_indios_sub$pop_ratio>=1 &agn_indios_sub$only_one==1)]))

strg  <- length(unique(agn_indios_sub$name[which(agn_indios_sub$per_encomienda>=enc_mean &agn_indios_sub$only_one==1)]))

weak <- length(unique(agn_indios_sub$name[which(agn_indios_sub$per_encomienda<enc_mean &agn_indios_sub$only_one==1)]))

dec_str <- length(unique(agn_indios_sub$name[which(agn_indios_sub$pop_ratio<1 & agn_indios_sub$per_encomienda>=enc_mean &agn_indios_sub$only_one==1)]))

dec_weak <- length(unique(agn_indios_sub$name[which(agn_indios_sub$pop_ratio<1 & agn_indios_sub$per_encomienda<enc_mean &agn_indios_sub$only_one==1)]))

unit_line <- c("Units", all, dec, inc, strg, weak, dec_str, dec_weak)

########
stargazer(model_dict2_all , model_dict2_pop_dec, model_dict2_pop_inc, model_dict2_le_strong, model_dict2_le_weak,
          model_dict2_pop_dec_strong , model_dict2_pop_dec_weak, 
          se   = list(model_dict2_all_vcov[, 2], model_dict2_pop_dec_vcov[, 2], 
                      model_dict2_pop_inc_vcov[, 2], model_dict2_le_strong_vcov[, 2],model_dict2_le_weak_vcov[,2] ,
                      model_dict2_pop_dec_strong_vcov[,2], model_dict2_pop_dec_weak_vcov[,2]
          ) ,
          p    = list(model_dict2_all_vcov[, 4], model_dict2_pop_dec_vcov[, 4], 
                      model_dict2_pop_inc_vcov[, 4], model_dict2_le_strong_vcov[, 4],model_dict2_le_weak_vcov[,4] ,
                      model_dict2_pop_dec_strong_vcov[,4], model_dict2_pop_dec_weak_vcov[,4]
          ),
          dep.var.labels= c("Favorable Ruling"),
          covariate.labels =c(  "Community", "Mistreatment", "Taxes"
          ),
          column.labels = c("All", "Population", "Local Elites", "Population Decrease"
          ),
          column.separate = c(1,2,2,2),
          object.names = F,
          #model.names = c("", "Decrease", "Increase", "Strong", "Weak", "Strong", "Weak"),
          keep = c("themeCommunity", "themeMistreatment", "themeTaxes"), 
          keep.stat = c("n", "adj.rsq"),
          add.lines=list(c("Province FE", "Y", "Y", "Y", "Y", "Y", "Y","Y"),
                         c("Time FE", "Y", "Y", "Y", "Y","Y", "Y","Y") ,
                         c("Province time trends", "Y", "Y", "Y", "Y","Y", "Y","Y") ,
                         unit_line
          ),
          label="only_one",
          title =("Favorable Ruling, Topics, and Balance of Powers (Single topic-cases)") 
#          out="only_one.tex"
)



#########################################
##  Supplementary material 
##  Acatlan and Piastla figure (Supplementary materials S7 )
###########################################

## Read datasey
acatlan <- read.csv("C://Users/franc/Dropbox/PhD/dissertation/paper_colonialism_liberalims/data/acatlan_piastla.csv")

p <- ggplot(acatlan, aes(x=decade, y=proportion)) + geom_line() + xlab("Decade") + ylab("Proportion of population under encomienda")
p <- p + geom_vline(xintercept = c(1550, 1570, 1600, 1650, 1745, 1800), linetype="dashed", color="gray50") 
p <- p +ylim(0,1) + theme_bw()
p
#ggsave("C:/Users/franc/Dropbox/PhD/dissertation/paper_colonialism_liberalims/tex/2019_version/plots/acatlan.png",p, width=6)

